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I. INTRODUCTION 


The problem of transmission of sound through shear layers has stimu- 
lated many investigations for quite some time. The simplest model con- 
siders a shear layer of zero thickness or an interface of two fluid media 
in relative motion known as the plane vortex sheet model. As early as 
1944, Landau^^^^ investigated the stability of such compressible fluid 
motions which evidenced a discontinuity in tangential velocity at a plane 
surface. He demonstrated the stability of the vortex sheet, that is, the 
persistence of infinitesimal displacements from its equilibruim position 
provided that there is a sufficiently large velocity jump across the dis- 
continuity. The nature of vortex sheet phenomena has acquired a greater 
importance in connection with theories of jet noise initiated by 
Lighthill^^^^ in 1951; specifically, it is appropriate to determine vortex 
sheet motions in the presence of an incident acoustical excitation. The 
problem investigated in this thesis concerns, in particular, the inter- 
action of an implusive ring source with a curved vortex sheet which forms 
the boundary of a cylindrical jet flow. 

During the late 1950 's numerous studies of reflection and trans- 
mission of sound by a vortex sheet due to incident time periodic plane 
acoustic waves were undertaken; after some initial efforts in which a 
kinematic boundary condition at the sheet was improperly formulated, 
Miles^^^^ (1957) removed this defect and successfully determined the plane 
wave reflection coefficient. He confirmed Landau's results regarding the 
stability of small amplitude vortex sheet motions and furthermore drew 
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attention to the so called neutrally stable motions with large amplitudes 
of the sheet that do not tend to grow or decay. The designation of 
resonant mode is given to such a large amplitude excitations of the vortex 
sheet. Miles' study was paralleled by that of Ribner^^®^ (1958) who also 
included transmission as well as reflection in a similar problem; both 
their studies were relatively simple and significant as prototypes of more 
complex problems. Miles^^^^ (1958) also analyzed the displacement of an 
initially undisturbed plane vortex sheet following a suddenly imposed, 
spatially periodic velocity disturbance. The resulting stability 
criterion was shown to be in agreement with that appropriate to the prior 
time periodic case. Transient wave propagation problems are inherently 
more difficult than the steady state harmonic ones and it was not until 
1969 that Friedland and Pierce^^^ determined the reflected field due to an 
implusive line source embedded in a stationary half-space, separated by a 
vortex sheet from a moving half-space. Howe^^^^ (1970) examined a similar 
problem, locating an impulsive point source close to the vortex sheet and 
stressed the fact that a correct analysis must account for both the insta- 
bility and resonance which are left out of the classical ray theory. 
Subsequently, in a series of papers published in the early 1970' s, Jones 
and Morgan^^^"^^^ and their associates^^^"^^^ investigated various 
problems of sound-vortex sheet interaction due to line or point sources 
with time-periodic or impulsive nature and significantly clarified the 
propagation of instability waves. 

Integral transform methods are employed for the analysis of the 
multi-variable (space-time) problems just mentioned; the technical diffi 
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culties» which arise predictably at the stage of inverse transformation of 
the solutions, are approached differently by various investigators. 
Causality can be ensured when the path of a time inversion integral is 
suitably located. This is accomplished for instance by Friedland and 
Pierce before analyzing the response integral with the help of Cagniard's 
technique which sequentially involves an integration variable change, a 
contour deformation, and an interchange of order of integration. Although 
they are able to derive an exact solution for the reflected wave, the 
physically more significant transmitted wave is set aside as being too 
difficult. Howe likewise adopts Cagniard's technique in his analysis and 
by restricting attention to the far field is able to deduce an explicit 
representation for the transmitted solution which incorporates the insta- 
bility and resonance waves. Jones and Morgan elected to invert the space 
Fourier transform first. The solutions thus obtained do not satisfy 
causality, that is, disturbances may precede the excitation, unless a 
homogeneous (source free) solution which corresponds to instability waves 
is added. This additional term is neither a conventional nor a real genera- 
lized function, but a delta function with complex arguments or so-called 
ultra-distribution. It is, in fact, the ultra-distribution thus intro- 
duced which helps to clarify the role of instability waves in a transient 
problem. However, the final form of their solutions is usually too 
abstract for extracting quantitative estimates of the excitations under 
consideration. 

More recently, Chao^^^^ (1977) employed the Cagniard's technique and 
also allowed for ultra-distribution in the inverse analysis; moreover he 


3 



was able to express in an exact manner the point source solution in terms 
of a single finite range integral and to characterize the solution without 
any integral in the line source case. The solutions which he arrived at 
are advantageous in that they enable the influence of a plane vortex sheet 
on an incident impulsive source of sound to be ascertained through uncomp- 
licated numerical evaluations. The outcome confirms the following. 
First, areas on the vortex sheet reached by the sound waves will now act as 
source regions which generate transmitted as well as reflected waves in the 
respective media. Second, for Mach number M < 2*^, sound disturbances 
will cause the vortex sheet to become unstable, resulting in the propaga- 
tion of instability waves. Finally, when M > 2, resonance or neutral 
stability waves are propagated. Both instability and resonance waves are 
shown to propagate downstream in a limited wedge-like region close to the 
vortex sheet. 

The success shown by Chao's method has motivated the investigation of 
a more complicated vortex/source geometry considered in this thesis. 
Specifically, the configuration chosen involves a curved vortex sheet 
which forms the boundary between an external medium at rest and a cylindri- 
cal jet column, otherwise known as an axisymmetric jet. And for mathemati- 
cal simplicity as well as practical interest, the source considered here 
assumes the shape of a circular ring (with negligible internal radius) 
whose plane is normal to the jet axis and is centered within the jet. Both 
fluids are assumed to have the same density and speed of sound to keep the 
details simple and all non-linearities including those of viscosity and 
heat conduction are neglected. 
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The stability of axisymmetric jets was first discussed by Batchelor 
and Gill^^^ (1962) and then in a similar fashion by Crow and Champagne^^®^ 
(1970), with criteria that involved not only the jet speed but also the 
wavelength of the primary sound source. Their works were significant in 
strongly suggesting that instability in jets is temporal, i.e. time-wise, 
in nature. The mathematical problem of cylindrical vortex layer insta- 
bility under the influence of a time periodic pressure fluctuation was 
first investigated byTam^^®^ (1971) but his solution simplified to that of 
a plane vortex sheet due to the nature of his approximation. However, his 
work was significant in linking the instability predicted by the vortex 
sheet model with the strong directional acoustic radiations from super- 
sonic jets observed in some experiments. Morgan^^®^ (1975) was actually 
the first to obtain a solution for the transmitted sound field in the 
motionless fluid outside an axisymmetric subsonic jet, due to a harmonic 
point source located off axis within the jet. He rigorously proved that 
there is essentially only one instability wave which arises, and that 
therefore an infinite number of other waves which are solutions to the 
homogeneous problem can be ignored. This investigation differs from that 
effected by Morgan in so far as a ring rather than a point source is used 
and furthermore that both subsonic and supersonic flows are contemplated. 

The organization and method chosen here for analysis are as follows. 
In Chapter II, the problem is formulated and its solution in the Fourier 
transformed space is obtained in the form of a multiple integral whose 
integrand involves complicated Bessel and Hankel functions. Chapter III 
contains the derivation of an approximate version of the integral solution 
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for the transmitted field; this is based on the assumption that the jet 
source radii are sufficiently large so that the two leading terms from the 
appropriate asymptotic (large argument) series expansions of Bessel and 
Hankel functions can be used. It is shown that proper attention is 
required on the different forms of Bessel function expansions over the 
argument regions. As a limiting case of both radii becoming arbitrarily 
large but retaining a finite difference between them, the first term is 
shown to reproduce the exact solution to the simplest model previously 
described, i.e. the transmitted field across a plane vortex sheet due to a 
line source embedded in the flow medium. The second term corresponds to 
the first order curvature effects. In approximating with these expansion 
terms, contributions arising from internal reflections are isolated and 
neglected, leaving a representation for only the primary or directly 
transmitted field. Howe^^^^ (1970) suggests that instability becomes sig- 
nificant after a time which is long enough for sound to travel backwards 
and forwards several times between opposite faces of the jet, when the 
acoustic coupling between these faces becomes large. The present analysis 
essentially calculates the transmitted field due to that sound which is 
emitted before these internal reflection effects become important. 

Singularities of the integrand are examined in Chapter IV prior to the 
inverse analysis. The branch points and the real pole which exists when 
M >2 associated with neutral stability wave are identical to that found in 
the plane vortex sheet case. However, the instability pole is shown to be 
a function of wavenumber as well as frequency in contrast to the plane case 
where it is independent of the former. Furthermore, as the pole remains 



complex for all Mach numbers, its transition to the neutral stability mode 
observed in the plane case does not take place. The coupling that exists 
between frequency and wavenumber in the integrand is responsible for this 
difference and indubitaly acts to complicate the inverse analysis. The 
condition of causality would be satisfied unambiguously if the double 
integral is carried out first with respect to frequency, but a procedure is 
not found to do this. Instead, in a fashion similar to Chao's application 
of Cagniard's technique, a change in the order of integration is facili- 
tated by deforming the frequency integration contour consistent with 
causality. In the plane case, this then leads to the simple solution forms 
alluded to previously. The solution obtained here is expectedly less 
simplex. Due to the curvature effects, the solution now exhibits a radial 
decay factor of 1//F~. In the evaluation for the specularly transmitted 
field, a branch point particular to the present problem is produced. As a 
result, in addition to the term that is similar to the corresponding field 
for the plane case, a second term arises from integration along the 
resulting branch cut. Evaluation of this branch cut integration is accomp- 
lished numerically by Gaussian interpolation. Related to neutral 
stability wave for M > 2 is a pole in the wavenumber space which is again 
distinctive to the curved vortex sheet problem as a consequence of its 
curvature being finite. This pole results in an additional term that is 
shown to decay exponentially away from the neutral stability wave and thus 
appears to effect a broadening of its singular wavefront. 

The integral that remains unresolved is one in which the moving insta- 
bility pole is directly involved. An analysis of this Integral is 
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presented in Chapter V where the essential features of instability wave is 
made relevant. First, it is shown that in the limit of kr^^ » where k 
is the wavenumber and r^ is the jet radius, the pole reduces to the 
instability pole for the plane case as anticipated. The pole is observed 
to move when kr„ is varied, but as kr^ becomes small, its movement cannot 
be adequately described by the approximate form derived from the two-term 
asymptotic expansions. Even though several terms would be needed to 
approximate the profile generated by the moving pole, it is shown that only 
one additional term is sufficient in establishing the general trend. A 
clear picture emerges when movement of the pole is determined for small kr^ 
as well by taking the leading terms in the convergent power series of the 
cylinder functions. From these considerations, the way in which the pole 
moves is approximated for all Mach numbers. Of particular interest is 
for M > 2 ^ where neutral stability instead of instability is predicted 
in the plane case because the complex pole becomes real. Here the pole is 
real only in the kr^ ® limit but remains complex for other values 
indicating the persistence of instability even for such high Mach numbers. 
As a consequence, M = 2*^ is not a clear cut transition point as 
predicted in the case of a plane vortex sheet that when M < the 
system is unstable and when M> 2/2~ the system is neutrally stable. In 
fact, instability waves may arise at all Mach numbers but further analysis 
shows that they are dependent on kr^ and cease to exist for sufficiently 
small kr^. In other words, instability waves will be caused by wavelengths 
large (high frequency) compared to the jet radius for all Mach numbers. 
The way in which this dependence takes effect is also shown to relate to 
the spatial extent of instability wave's presence, which is seen to con- 
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tract towards the vortex sheet with increasing wavelength. This predic- 
tion appears to be in agreement with the strong directional acoustic radia- 
tion observed in some jet noise experiments. In Chapter VI, some numerical 
results are presented and discussed. 

The idealized problem of the curved vortex sheet enclosing a 
cylindrical column of jet presented in this thesis is of particular 
interest since it is a first approximation to the real jet. As a next step 
in the modelling, it is desirable to include a rigid attachment to the flow 
to simulate a jet nozzle. This will bring out any contribution to the 
noise generation process of edge scattering which may in fact be extremely 
important in the overall picture. Problems associated with edge 
scattering have mixed boundary values and are in general rather difficult 
to handle. Nevertheless, a significant literature on that subject exists, 
including contributions by Crighton^^^ (1972), Crighton & Leppington^^^ 
(1974), and more recently by Miura^^^^ for the plane case and by Munt^^^^ 
(1977) for the curved case, just to mention a few. It is also desirable to 
consider a finite shear layer instead of an infinitesimally thin vortex 
sheet since in all practicality it never occurs even near the jet nozzle 
exit because of boundary layer. Some investigators have ventured into 
plane shear layer models such as Graham & Graham^^®^ (1968), Jones^^^^ 
(1977), and Koutsoyannis^^®^ (1980) for example, but none has considered a 
curved shear layer. 
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II. FORMULATION OF THE PROBLEM 


Consider the problem in which a cylindrical column of jet fluid moves 
axially through a similar stationary fluid medium, extending over the 
entire space. The two fluids are ideally separated by a vortex layer that 
is equivalent to a shear layer of zero thickness. The jet is taken 
with its axis along the z-axis of a cylindrical polar coordinate system 
(r, 6 , z) and to have constant speed U and radius r^. Within the jet, 
oriented along the z 0 plane, is placed a ring shaped mass source of radius 
Tg as shown in Figure II-l. 


r 



FIGURE II-l Sketch of the Problem Geometry 

Without affecting the underlying phenomena, both fluids are assumed 
to have the same speed of sound and density to keep the details simple. 
All nonlinear effects including those of viscosity and thermal conducti- 
vity will be ignored, which is expected to introduce some uncertainties, 
especially in the amplitude prediction of instability waves. Because of 
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the ax 1 symmetric nature of the problem, all quantities are independent of 
0 variation and the equations satisfied by the velocity potential 
4> (r, z, t) are then as follows: 


in region I outside the jet 




3r 


- = 0 


r > r. 


(II-l) 


in region II inside the jet 


■M-rr ^ “ -fe-) ^*2 - = 


Q«(2)S(r - r )6 (t) 

^ r< r 

2 nr_ ° 


where 


V 





(II-3) 


a is the speed of sound and Q is the strength of the ring source which is 
impulsively activated, i.e. turned on and off at time t = 0. The plane 
vortex sheet geometry is a special case of the above when r^ and r^ become 
arbitrarily large. Then, by the fact that Q = 2 7rr^q, where q is the 
source strength per unit length. Equations II-l & 2 are replaced by: 


a'^ at ^ 

y >0 

(Il-la) 

V(n ^ “ L)'*2 -’'*2 

d 

y < 0 

{II-2a) 

■ q 5 (x) d (y + y^) 6 (t) 
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where now 


7 


2 





+ 



(II-3a) 


X being the axial coordinate of the rectangular coordinate system (x,y,z), 
the vortex sheet being located at y = 0, and the line source at 
y = -y^, x = 0. This line source problem had been investigated by Chao^^^ 
and his results will form a basis for comparison with the present problem. 


There are two matching conditions at the common boundary r = r^ of the 
two regions of irrotational motion. The dynamical condition that the 
pressure is continuous across the boundary gives 


8 4' <1 9 '3<f> o 

+ U — 

3t 3t ^ 3z 


(II-4) 


The kinematic condition for equal particle displacement n on both sides 
of the vortex layer gives 


_ 5 n 
3r " 3t 


^'^2 _ 3 H . AH 
TT ^ 3 z 


(II-5) 


In addition, there are the following conditions that must be 
satisfied by the solution. The radiation condition requires that the field 
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be radiating outwards at infinity, the finiteness condition is required at 
the center of the jet that insures no blow up of the solution there, and 
finally the causality condition, which states that the field be unper- 
turbed before time t = 0, will apply. 


Define (r,k,u»), the Fourier transform of 

;.7: .j • " “ 

's' LL « *• * 


4> (r,z,t) as: 


(II-6a) 

(II-6b) 


where j = 1, 2, and Cj^ and are integration contours in the wave number 
k-plane and the frequency aj-plane respectively which will satisfy the 
causality condition. 


Then Eqs. II-l & 2 reduce to 


( 

( 



1 d 
rdr 


)’fi 


j. 

r dry 




2 

+ Yj 


1 


+ Yj 


0 


Q5 (r-r^) 

o 3 ^ 

8 Ti r^ 


where 



Y 




(II-7) 


(II-8) 


(II-9) 
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M = U/A is the Mach number of the flow. 


The boundary conditions at r = r^, Eqs. II-4 & 5 become: 


(1) = (u - Uk) ^2 


(II-IO) 


(w- Uk) 


d ipj 

dr 


(0 


d il> 2 
dr 


(II-ll) 


Eqs. II-7 & 8 can be recognized as the zeroth order Bessel's equations 
and the solution to Eq. II-7 which is for region I is 

'*'1 " ^1 (Yj r) + Bj ( Yj r) (11-12) 

where and are Hankel functions of the first and second kind 

respectively, and B^ are constants to be determined. The y's are 
generally complex and if the argument of Y is restricted to non-negative 
values, i.e., Im(Y > 0, then by the radiation condition, B^ = 0 and 
^ 2 »*educes to 

' '' i ( ■^1 '•) 




(11-13) 


For region II, the solution to Eq. II-8 can be put down as: 



(11-14) 


where and are Bessel's functions of the first and second kind respec- 
tively, and as the equation has an inhomogeneous right hand side, A2 and B2 
are functions of r determined by the variation of parameters. Thus, 



Q 

1671 ^ 

_ _Q 

1677 ^ 


Yo 2 ^ - »*s) + Cj 

Jo (Y 2 H (r - r^) + C2 


(11-15) 


where and C2 are constants to be determined and H stands for Heaviside 
function. 


t 

It is to bo notod thdt in tho limiting caso of Y ^^00 

However, with H(r - r^) = 1 - H(r^ - r), Aj can be rewritten as 

''2 ' • 2 '‘o 2 ''s' " <’■5 - >•) (II-16) 

where now as r. h- O, A„ c! . 

* t 1 


X5. 


By the finiteness condition, C£ is necessarily zero since as r 0, 


Yp ■+■ — 00 , Therefore, Eq. 11-14 becomes 


h " ^1 *^0 ^ ^2 2 


(11-17) 


where S^(r) = ( Y 2 r^) 0^ ( Yg r) H (r^ " »*) + Jq ^^2 

Even though the condition Im( Y 2 ) >0 is not required, that branch is 

assumed for consistency with Y^ and this choice also happens to yield a 
better form. The constants Aj^ and Cj can now be solved by applying 
boundary conditions Eqs. II - 10 & 11. 

Q (o( (0- UK) Jp (Y 2 r^) 

Btt^ r S, 

0 J 


g ■Jq ■'s) % 

I 61 T ^ Si 


(11-18) 


where Sj = Yj( UK)" ujH (K^r^) (Y^r^) OjCf^r^) 

Sy UK)" Hp> CTjrJ V„ (r^r^) - Cr_,r„) Yjjfr^r^) 


Finally, the transformed solutions are obtained as: 


'i>^ = 


Q (u((o- UK) Jp (Yg r^) (Y^r) 

8tt ^ r„ S, 
o 0 


(11-19) 



and 


= 


16 Tr‘ 




Jo ("z ’■o) Jo '•) Sy 


- s„ (-•) 




(11-20) 


The purpose of this work is to derive for the transmitted field 
exterior of the jet alone, and therefore only Eq. 11-19 will be dealt with 
for its inverse. This, in reality, is clearly of more practical value and 
interest than the corresponding reflected field internal to the jet. How- 
ever, one may seek to find it as well by inverting Eq. 11-20 in a manner 
similar to that which will be described in the following text, although 
more difficulties are expected, as if; 2 appears substantially more complex 
than tpj. The integral in question is given by Eq. II-6b and substituting 
Eq. 11-19 into it gives: 


= _ 5 _ 




-ff 


»+iE io(a.-UK) 


0 


( 11 - 21 ) 


A brief discussion on the selection of integration contours is appro- 
priate here. Note that integration with respect to time is to be carried 
out prior to space. This is to ascertain that the causality condition can 
be dealt with rigorously and unambiguously. If space integration was to be 
performed first, then the contour Cj^ does not have a clear cut choice and 
one may obtain a solution which is non-causal. The way in which causality 
is satisfied is by selecting an c large enough such that any and all 
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singularities associated with the integrand will lie below the contour. 
This requires that none of the zeros of Sj be at infinity, an assumption 
which will be shown to hold in a later section. It is then apparent that 
the field is quiet before time t = 0 by invoking Jordan's Lemma. Since the 
causality condition is duly fulfilled in the time domain and there isn't a 
comparative condition in the space domain, Cj^ assumes the usual intergra- 
tion contour along the real axis. 

0 

The proposed method of attack is in the spirit of so-called Cagniard's 
method which was in wide application on seismic wave propagation problems 
until recent times. The method may be summarized in simple terms as a 
sequence of a change of integration variable, a contour deformation, and an 
exchange in the order of integration. Friedland and Pierce^^^ were first 
to apply this method on any problem involving wave propagation in fluid 
mediums. Its success depends largely on the particular contour chosen for 
deformation, as exemplified by Chao^^^. 


First, make the change of variable. Let 
4 = ^ (11-22) 


which is the complex phase speed. This differs from the original 
Cagniard's method in which k/u or the wave slowness is the new variable. 
Eq. 11-21 becomes 
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£(C-M)Jo(l(rgr^) H^^)(kr,r) 




k(^at-z) dk 


(11-23) 


where 

V -r2E^H<')(krjrj,)Jj(kr2r^) (11-24) 

and 

kr^ = k (^2 _ jjJs = 

kr2 = k[(c-Mj^ = 

The contour C ^ = C(k) is now a function of k in accordance with Eq. 
11-22 and is shown below. It consists of two horizontal paths such that all 
the singularities lie below when k > 0 and above C_ when k < 0. The 

causality condition is observed by again applying Jordan's Lemma. 



^ - plane 



Figure II-2 Integration Contour C(k) and 
Branch Cuts for r^, Tg 


r 2 and T 2 defined by Eq. 11-25 have four branch points at 
C=+l and C=M±1 

To satisfy the commitment made previously which restricts Im( Yj) > 0 and 
Inj( T 2 ) ^ 0, the imaginary parts of r ^ and T 2 must have the same sign 
as k while integrating along C(k) in the ^ -plane. Accordingly, branch 
cuts are placed as illustrated in Fig. II-2. This restricts each of the 
arguments of C ± 1 snd C i 1 ^ between the range of - ir and while 

the range of argument for both and F 2 will faill within 0 and n 

above or -tt and 0 below the real axis for all points, thus automatically 
satisfying the requirement on the F 's. Due to the complex nature of Eq. 



11-23, an exact expression for 4) is not possible and some approximations 
will have to be made next. 
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III. THE APPROXIMATE SOLlfTION 


a. ASYMPTOMATIC EXPANSION OF BESSEL FUNCTIONS 


Carrying out the exact inverse transform as presented by Eq. 11-23 has 
not been possible because of the complicated Bessel and Hankel functions in 
the integrand. These cylinder functions are now approximated. First, 
observe their arguments to contain r^, the jet radius, r^, the source 
radius, and r, the radial distance out from the jet. These are physical 
parameters which can and will be assumed to be large for the present 
problem. Then, provided kr^ and k T 2 remain finite, the asymptotic (large 
argument) series expansions of the cylinder functions can be used, r ^ 
and r 2 can never be zero since all their zeros must lie within the strip 
bounded by the integration contours C^. But the integral in k clearly will 
pass through zero and the arguments will become vanishingly small near that 
point. However, it will be shown that contribution from small k can be 
ignored except in dealing with instability. The expansions for Hankel 
functions are: 



(III-l) 


But the expansions for Bessel functions take on different forms. 
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where plus is for — tt < argz < ir and minus for 0 < arg 2 < 2 . 

While asymptotic expansions of functions with complex arguments may vary 
in different regions, they usually do not take on different forms in an 
overlapping region as indicated by Eq. III-2. That the expansions for 
and differ by the small term in the range of argument between 0 and ^ 
is recognized as Stokes' phenomenon, so called because it was first dis- 
covered by Stokes^^^^ in 1857. It has been confirmed in the present study 
that a correct analysis must account for this subtle difference. This 
situation does not arise if only one term in the expansion series is used 
as Tam^^®^ did in his analysis. But the solution then simplifies to that 
for the plane case because curvature effects enter through higher order 
terms. The two forms of expansions for J (z) are each applied in the 
following manner. Recall in Chapter II the branch cuts of and P 2 
are defined so that their arguments fall within - tt and tt. This means 
for k > 0, 

-7T < arg (kr.) < tt 

but for k < 0, 

0 > arg (kFj) <; 2ir , j = 1, 2 
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By approximately combining Eq. III-l with Eq. III-2, the following 
approximation is obtained. 


Oo (kr^Fj) H<l)(krFi) = !■, C“ 


s'r2 




^^ik(rrj-r^r2).E^,1k(rrj+r^F2)]1 

.1kro(rj-F2) I I 


Jj (kr^r^) H<^>(kr^Fj) = 






(III-3) 


where plus sign is for k < 0 and negative sign for k > 0. 


In the above. 


Dj = 1 + 


l—( _i L ) 

8k V r^r2 rfj ^ 


8k 


V“( "f~)“ ■' 

•^0 ^2 ^1 






Tz h 

1- (-J- + 

8k 


ihy^ 


Eo = 1 


. 4 (_i ^ ) 

8 R ?7 ^2 Tj ^ 


(III-4) 
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Now, to show that contributions from small k can be ignored, leading terms 
from the power series expansions of the cylinder functions are substituted 
into the integrand of Eq. 11-23. 


kr r, 

E(e-M) 

0 1 *• 


) 




ln( 


kr. 
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(III-5) 


In the limit of k -i- 0, the above quantity can be seen to vanish except at 

the zeros of the denominator. When treating with instability, its effect 

on the integral will be studied. Here, Eq. 11-23 is replaced by the 

approximate form using Eqs. III-3 and 4. 
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(III-6) 

d^dk 


where 


f = 2 + Tj (r - r„) + T2 (r„ - r^) 5at 


(III-7) 
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i ,22 (III-8) 

^2 = (^1 -^z) + srjFg) - (Fj - Tg) (1 -rjFg) 

This expression represents the transmitted field from a large ring-shaped 
source positioned within a cylindrical jet of large radius r^ with the 
source plane normal to the jet axis. Hence, the radial field distance r is 
at least as large as r^. In the special case of r^ and becoming 
arbitrarily large but retaining a finite difference between them, Eq. III- 
6 reduces to 



(III-9) 

which is not unexpectedly the exact solution to the line source problem as 
stated by Eqs. Il-la through 3a. Hence, it may be concluded that the small 
k or high frequency effects may be lost when a real jet is approximated by 
a plane vortex sheet. 

b. FURTHER APPROXIMATIONS 


The approximate expression of Eq. III-6 for the transmitted field 
4> 2 is still too complex for evaluation and further approximations are 
necessary. An examination of the integrand shows its denominator to con 


tain a term exponentially small due to the fact that ImtkTg) > 0, and 
this permits the following series expansions: 


L Oj J Oj 

KTl) ' “ ^ ^ •••• 


(III-IO) 


Rewrite Eq. III-6 as 








I f°r r,^ gi'2krjr2 ^ ^ e'^^'^0^2 +....T 

\l i- °1 “1 “1 


g ’2kr,r2 ,12kr^f2 . ,],1kf 


(III-ll) 


Thus, the integrands in Eq. 11-23 have been transformed to infinite series 
following two successive approximations. The integrals are assumed to be 
uniformly convergent so that the order of integration and summation can be 
interchanged. Eq. III-ll then represents the transmitted solution as a sum 
of infinite number of integrals. By attaching physical significance to 
contribution from each integral, it will now be demonstrated that only a 
finite number of terms are needed to generate a fairly descriptive picture 
of the transmitted field. First, it is simple to observe the exact solu- 
tion to the line source problem as given by Eq. III-9 is again reducible 
from Eq. III-ll, requiring however only its first term. Second, a wealth 
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of information can be gleaned through the examination of exponential fac- 
tors in these integral solutions. For instance, this factor in Eq. III-9 
carries with it the sense of a wave front and subsequent motion initiating 
from X = 0, y = -y^, which are points on the line source. This reasoning 
will be substantiated in a later section for the factors found in each term 
of Eq. III-ll, which have a parallel functionality. The first of these is 
f as defined by Eq. III-8a and it is identical to the corresponding factor 
in Eq. III-9 just described. Fig. Ill-la illustrates the fact that f 
carries the information on wave front as it relates to signals coming from 
the source point closest to the observer field position. The factors in 
the subsequent terms all contain f and an expression based on k T2» and it 
is not difficult to find an interpretation for each of them. Thus, the 
second term relates to signals coming from the source point farthest to the 
field position as depicted by Fig. Ill-lb. Since the signals will be 
continuous following the wavefront, combination of the first two terms of 
Eq. III-ll is seen as descriptive of the transmitted field due to the 
primary source-vortex layer interaction, beginning with the wave front 
coming from A and ending with the contribution from B. The primary inter- 
action is clearly without any reflected contributions, which as expected 
can be found in the succeeding terms. The field resulting from a single 
reflection is contained in the third and fourth terms as indicated by their 
respective exponential factor (fourth term is not shown in Eq. III-ll but 
basically it is the product of second and third terms). As illustrated in 
Figs. III-lc and Id, it begins from A and ends at B with a single reflec- 
tion at the interface. This simple argument can be applied on subsequent 
terms to show that they represent multiple reflections of increasing 
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order. suggests that after several reflections, the acoustic 

coupling between opposite faces of the jet may become large. It 

X Any Field Point x x x 

i ^ ‘ 

w® ® 

a. 1st term b. 2nd term c. 3rd term d. 4th term 

Figure III-l Physical Significance of Terms in Eq. III-ll 

will be assumed here that for small time the contributions arising from 
these internal reflections are negligible to the total transmitted field. 
With these considerations, just the first two pair of terms in Eq. III-ll 
will be evaluated to focus attention on the added role played by the 
curvature effects in the transmission mechanism. 

c. NON-DIMENSIONAL FORM 


First, rewrite Eq. III-ll as follows. 
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where 


g = z + Tj (r - r^) + T 2 (r^ + r^) - 5at 


Notice the similarity as well as the difference between and <l> , They 

9 

will be treated separately from here on. 



d^dk 


(III-13) 


where 


Aj = (Tj + (1 + Fj Tg) 


(III-13a) 


Aj = i - r2^)(l - 3 Fj r2)/8 
1 - iro / 1 1 X 

--8 ^ T2r^ -^7-T^J 


(III-13b) 


A check on Eq. III-13 reveals the imaginary part to vanish as anticipated 
and hence it can be put in the following equivalent form, with Re to mean 
the real part of the integrals. 

aQ Re 
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(III-14) 


After some algebraic manipulations, this equation is non-dimensional ized 
as follows: 
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(III-15) 


with non-dimensional quantities: 


P »R“^ »K = kr^j » Z = ^ 


0 'I 
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<f>g is similarly put in its non-dimensional form: 
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Here, 
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Z + Fj (R - 1) + Fg (1 + R5) - Ct 
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^*P and do not differ greatly in their appearance but retracing from 
previous discussions on the significance of each term, it is to be kept in 
mind that while the former contains information on the wave front, the 
latter should be free of it since there can only be a single wave front. 
Whether or not this is indeed the case will become apparent in the final 
evaluation of the solutions. 
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IV. INVERSE ANALYSIS 
a. THE CONDITION OF CAUSALITY 


In two previous chapters, an integral solution of the transmitted 
field which approximates for large values of radial lengths and which 
excludes any contributions arising from internal reflections is derived. 
This inversion integral with its non-dimensional form as given by Eqs. III- 
15 & 16 is evaluated in the present chapter to the simplest analytic repre- 
sentation possible. To begin with, an examination of the integrand is made 
(see Appendix A) to identify all the singularities. The branch points and 
the real pole which exists when M > 2 associated with neutral stability 
wave are identical to that found in the plane vortex sheet case. The 
instability pole is complicated by being a function of wavenumber as v/ell 
as frequency in contrast to the plane case and detailed examination will be 
presented later in its appropriate place. 

Recall the integration paths for 5 and K are originally as follows. 
K runs from 0 to ~ while K moves along the path C^ such that all 
singularities are restricted to lie below it. Then, the condition of 
causality can be checked by invoking Jordan's Lemma as illustrated in 
Figure IV-2. 
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Since there is no singularities within the closed contour of C+ and 
the deformed path C , the integrals vanish by Cauchy's Integral 
Theorem. That is. 


R - Rg > T, = 0 

R + Rg > T, = 0 


^ (IV-3) 


The first implies that before time t equals R - R^, which is the minimum 
time for the waves to reach the field point R, $ p remains zero to satisfy 
the physics that before the secondary sources on the interface are excited, 
there cannot be any response from the transmitted field. Before elabo- 
rating on the second condition given in Eq. IV-3, it is expedient to 
reflect on the differences between <I>p and that are easily discernible 
in Eqs. III-15 & 16. The development of Chapter III pointed to the fact 
that their sum characterizes only the directly transmitted field, having 
disregarded all successive terms arising from one or more reflections at 
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the interface. A question comes up as to why the characterization goes in 
pairs, whose answer may very well rest with the curvature variations. For 
if the field point of interest is at P, as illustrated in Fig. IV-3, the 
top half of the circular interface is seen to have the curvature, say, in a 
positive sense. Top half of the ring source similarly has its curvature in 
the same sense, but that of the lower half clearly has the opposite sense. 
This makes it plausible for an interpretation of and in the manner 
described below. is representative of the transmitted field arising 
from top half of the source as depicted in Fig. IV-3b, and the first 
condition of Eq. IV-3 which implies the initiation of the wave front from 
point A nearest point P, applies as it gives the minimum time necessary for 
the waves to make contact with the field point. Obviously, transmission 
cannot occur prior to this time. On the other hand, 4' g represents the 
field originating from lower half of the source as is apparent from Fig. 
IV-3c, and the second condition of Eq. IV-4 which cannot have anything to 
do with the wave front, applies as again it furnishes the minimum time 
before transmission from any point on the lower half source is possible. 



Figure IV-3 
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It can be concluded from the previous discussion that the integral 
solutions represented by Eqs. III-15 & 16 exhibit the causality conditions 
properly and this clears the way for the inverse analysis. For the 
inequalities of Eq. IV-3 to hold in the opposite sense, i.e. < instead of 
Jordan's Lemma requires taking the lower semi-circle on Fig. IV-2 which 
now together with C+ forms a region containing all the singularities as 
well as the branch cuts. Basically, the inverse analysis entails the 
evaluation of residues coming from the singularities and integrating 
around the branch cuts. While residues are relatively simple to evaluate, 
branch cut integrations for complicated integrand such as the ones 
encountered presently are not likely to be easy nor manageable to yield a 
final expression useful for applications. This is the reason why a partic- 
ular technique is demanded for the analysis. 

b. ZERO IMAGINARY CONTOURS 


Cagniard's technique will be employed in similar fashion as Chao 
applied it to his plane vortex sheet problems. The success to this method 
lies in finding a path for contour deformation such that integration along 
the branch cuts need not be concerned with and an exchange in the order of 
integration can be made consistent with the causality condition. In other 
words, the idea is to find a way to integrate with respect to K first 
instead of C because that is substantially less cumbersome to carry out 
than the other way around. An examination of Eqs. III-15 & 16 shows that 
if a contour on which Im(F) or Im(G) is greater than zero can be found, 
then the power raised to the exponent in either integral will have a 
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negative real quantity for all K, and along such a path, the order of 
integration can be altered without any violation. So the functions F and G 
will be studied next. Recall 


F = Z + (R - 1) + Tg (1 - Rg) - Ct 

G = Z + Tj (R - 1) + Tg (1 - R 5 ) - 5 t 


(IV-4) 


First, it is observed that they are purely real along the real axis not 
covered by the branch cuts of and r^. A parametric equation which 
defines the curves of Im(F) = 0 or Im(G) = 0 not on the real axis can be 
determined. Chao showed in his plane vortex sheet problem that these 
curves are oval-shaped. Here the same is anticipated as the functions are 
quite similar, but the curves which shall be designated CpQ and ^GO are 
numerically searched instead. A computer program can be set up with 
relative ease and generate a whole array of them by varying the parameters 
X, R, Rg, and t. Some examples are shown in Fig. IV-4 with just the top 
half as they are symmetrical about the real axis. Several observations can 
be made regarding these curves. Whether or not they split into two depends 
primarily on the Mach number. For M < 2 when the branch cuts overlap, a 
single curve is maintained, whereas for M > 2 after the branch cuts become 
separated, two similar curves will evolve around the cuts when the time t 
is at least greater than that shown in the inequality of Eq. IV-3. When 
time is small, the oval-shaped curve is very large for all Mach numbers. 
As time grows larger, the curve starts to move in and becomes 
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FIGURE IV-4 EXAMPLES 0F ImF-0(L.) AND ImG-OCR.) FER Rg-O.! 
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smaller until when t gets infinitely large, the curve(s) will just barely 
wrap around the branch cuts. Finally, the curves usually bias to one side 
or the other depending on the values of R and as they resemble weights 
to and T 2 respectively. All the background information just secured 
on CpQ and Cqq will help in the analysis later. 

Now that the question of where and how the imaginary parts of F and 6 
vanish has been thoroughly resolved, it is possible to establish the con- 
tours along which the individual imaginary parts are always greater than 
zero. These shall be designated by Cp^ and respectively and a sketch 
is shown in Fig. IV-5. 



Figure IV-5 Sketch of Contours ^p+> 
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The plus and the minus syiibols in double prime quotes are indicative of the 
signs of F or G in the various regions. The jumps across contours where 
the imaginary parts vanish are distinctly clear. Cp^, or identically Cg^, 
for that matter, wades through the positive imaginary regions in the manner 
shown. Fig. IV-5a illustrates the case when M < 2 or M > 2 but t small. 
In such a case, the contour is inside of CpQ above the real axis and 
outside of CpQ below. When M > 2 and there are two of the zero imaginary 
curves as illustrated by Fig. IV-5b, then the contour takes a little detour 
into the lower half spaces between the branch cuts such that the local 
portion of the real axis outside of Cpg is exterior of Cp^ as well. A point 
which is rather self-evident will be made here. Note that whereas the 
contour Cpg or Cgg is uniquely determined for a given set of parametric 
values, there is no clear-cut choice on Cp^ or Cg^. This is simply due to 
the fact that regions of positive imaginary part spread over a large area. 
The added flexibility thus acquired in the selection of a contour might be 
interpreted on first thought to give rise to some technical difficulties in 
performing the actual contour deformation. However, taking a closer look, 
one would expect the final solution to be independent of the particular 
contour selected, because there cannot possibly be more than one solution 
for a problem so well defined as the present. This supposition will in 
fact be shown to be true. 

From here on, the analysis will be described for 4>p only, with the 
understanding that results for ^ g could be similarly obtained. It is easy 
to see that while some singularities will be enclosed by Cp^, others will 
be left out depending on the particular contour chosen. The final results 





however will be the same regardless of the choice made. Actually, the 
question of where the poles are located with respect to CpQ is of greater 
importance and is examined next. As the imaginary part of F is zero along 
CpQ, it is natural to look at the sign of F to see whether a point is in or 
out of the closed contour. First, consider q and ^ g which are given 
in Eq. A-2. Substitute them into F defined by Eq. IV-4 and sort out the 
real and imaginary parts. 

F = Z + + 1 - l] (R + Rj - 2) - 

± [(m2 (^.1) -t] 

= i B (IV-5) 


In agreement with Fig. IV-5, 2 will be inside of CpQ if 3 > 0 and 

outside if B < o. Using the Heaviside step function, let 
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1 



(R - Rs) 
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( Vm^+i'+ 



(IV-5a) 


E 3 » 54 ? 5 defined in Eq. A-2 are on the real axis along which 
the imaginary part of F vanishes. Therefore, some other check is needed to 
determine where they are in relation to CpQ. First of all, these poles 
will be enclosed by <^FO as long as it is a single curve. When they split 
into two, then a condition exists for the poles to fall outside of CpQ. 
Since its derivation is similar to that found in Appendix 1 of Chao^^^, 
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only the condition itself will be presented here. The poles ^345 
will be outside of CpQ if the inequality represented by the following is 
satisfied. 


V 



T > 




+ 


K - - '’s) 


m = 3, 4, 5 


(IV- 6 ) 


which is the pole of K + A 2 /( f j f 2 ^ 1 ^* examined 

later for its movement as a function of K and how that will relate to CpQ. 

c. EVALUATION OF INTEGRALS 


To begin with, break up 4> p given by Eq. III-15 into three parts as 
follows: , 
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(IV-9) 


(IV-10) 


$ g and $Q can be combined but they are evaluated separately for mathe- 
matical simplicity. As discussed at some length earlier, the main idea in 
the inverse analysis is to deform from to Cp^ as illustrated in Fig. 
IV-7. By Cauchy's Integral Theorem, integration along is now equivalent 
to the sum of integrations along Cp^, as -> ® and Cj's which 

gives the residues of all the poles exterior to Cp^. Because the 

2 

integrands decay as good or better than 1/ C for ^ becoming very large, 
the contour makes no contribution by Jordan's Lemma and thus Eqs. IV-8 to 
10 following the contour deformation becomes (with the poles indexed in 
accordance with Appendix A): 
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Figure IV-7 Integration Contours 
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In the representations above, an exchange of order of integration has been 
effected along Cp^ where the real part of iKF is positive. For the 
present, all poles have been assumed to fall outside of Cp^ even though the 
final outcome will not depend on where the poles lie in relation to Cp^. 

First, let 




iKF 


dK 


K + 


Fir^Ai 


(IV-14) 


and evaluate this integral which occurs in and ^ By a simple 
translation and a rotation of the K-coordinate, it is possible to transform 
Ij^ into a form involving exponential integral as follows: 


\ E. (- i K^ F) 


(IV-15) 


where K^ = - A 2 /( f i ^2 Aj) is the 
discussed in Appendix A and (-iK ^ F) 
integral 


moving pole in the K-plane 
represents the exponential 
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Thus, 
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Notice that $ and <i> are the same integrals but opposite in sign to 
drop out from the above. Now, the integrand of behaves as 



Fig. IV-8 Contour Deformation of 
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Integration over does not contribute as can be seen from the above by 
taking infinitely large. Thus, 



(IV-18) 


where Cq is evaluated if and when the zero of F lies exterior of Cp^. Cj's 
prevail by the previous assumption that the enclosed poles are all outside 
of Cp,. By the residue theorem, Eq. IV-18 is evaluated to give 







S'") Qp (5„) 


Co (Sq-H) 
^1 


The zero of F is denoted by ^ q and deserves some discussion for its 
physical significance. Consider first the fact that those points where F 
vanishes are also the intersection points for Re(F) and Im(F) where they 
both vanish. Earlier, it was shown (Fig. IV-4) that curves of Im(F) = 0 
which are denoted by CpQ are oval-shaped. In Fig. IV-9, the curves of 
Re(F) = 0 are illustrated along with the corresponding curves CpQ. These 
curves in general must and can be determined numerically without much 
difficulties. Basically, the real part of F is shown to vanish along a 
curve that approaches a straight line perpendicular to the real axis when 
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far away from ho- Only the upper half is shown in both Figs. IV-4 and 9 
because the curves are symmetric about the real axis. Furthermore, Cp^ is 
such that it always encloses the lower half of CpQ together with the branch 
cuts and hence the need to include them in the analysis is circumvented as 
a consequence of the contour deformation. 

The position of is found to be a function of time. For small 
time, it is on the real axis exterior of CpQ which when substituted into 
Eq. IV-19 gives rise to a purely imaginary quantity. This clearly does not 
contribute to the real transmitted field represented by ^ p of Eq. III-15. 
X Q moves toward CpQ with growing time (Fig. IV-lOa) and reaches the point 
where CpQ meets the real axis (Fig. IV-lOb). Beyond this time, ^ q 
becomes complex by moving onto and oscillating back and forth on Cpo (Fig* 
IV-lOc) for all times. This is true with one exception that occurs for 
M > 2 when CpQ splits into two parts as exemplified in Figure IV-9. For 
these cases, 5 q may fall onto the real axis between the two parts so that 
even for time greater than ^Qp, it can become real again. Whereas K q 
does not contribute for time x< iQp because it gives rise to a purely 
imaginary quantity, it always yields a real quantity to Eq. IV-19 for 
T> Tqp and therefore contributes to the transmitted sound field. Recall 
the condition of causality given by Eq. IV-3 which states that tp = 0 for 
T< R - R^, the minimum time of sound propagation from source to vortex 
layer. Through numerical calculations, is found to be the minimum time 
of wave arrival at a field point of interest and Tqf R - Rs obvious. 
Thus, provides just the geometric acoustics contributions. In fact, 
it has been verified that all the remaining contributions, which arise from 
sound- vortex sheet interations, arrive at later times than 
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Figure IV-10 Cq Motion Analogy of Wave Propagation 


Next, focus attention on ^>^2 which can be rewritten from Eq. IV-11 
as 
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When j = 1 and is inside of CpQ so that Im(F) > 0 on Cj^, 
in the order of integration is possible. Then, by using Eqs. 
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where for general purposes 
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For the case of j = 2, inside of Cp^ when its complex 

conjugate lies within Cpg. Thus, $ ^22 " 0 = 1 but is 

outside of Cpg so that Im(F) < 0 on Cp $^21 evaluated in the 

following manner. 
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(IV-21a) 


K + 




Eq. IV-21 alone is sufficient to represent 4’ for both cases since the 
Heaviside function which is defined by Eq. IV-5a, is zero when C ^ 
lies outside of Cpg and Eq. IV-20 is recovered. When is outside of 
Cpg, so is ^2 lie on the outside of Cp^. Then, 
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and ^ g are points on the real axis and exterior of Cp^ so that the 
imaginary part of F is greater than zero on the portion of and below 
the real axis and less than zero above. Thus, 
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In ^ A2Km’ integral along is for the upper half only. When j = 3, 
the situation is similar to j = 4 or 5 because ^ 3 is also a point on the 
real axis exterior of Cp^. 
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Evaluation of ***^26 saved for a later section as it requires 

special treatment. Note the absence in Eq. IV-25 of any additional factors 
to indicate when and where each term is present. For example, is 

present only when M < 2 V~2~and where it is exterior of Cp^ so that the 
corresponding terms would vanish if these conditions are not met. The same 
in a consistent manner is implied on Eq. IV-19 but the factors are not 
tagged on in the meantime for simplicity. By duplicating the procedure of 
^ following results are obtained. 
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Again, ^ will be evaluated later along with ^ [^ 2 ^' Next, 4>/^2K3 
<I> g 2 K 3 3 re simplified further by directly substituting ^3 “ M /2 into 
these integrals and working out the double integration. 
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= 2tt' 
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Qa(? 3 ) = ^3(^3-^) K3 




^A^^3^ it 


^1^2 


3M ? 8 


(IV-28b) 




^ 2^2.4' (8- m 2) 

^KF " ^[''(^3^] ^ [^"2] M [ 21 / 2 ”- m] - H [-F(C3)] H [m -2^^.] 

It will now be shown that ^ /\2Kj's* j = 1» 2, 4, 5 all vanish in 

reference to earlier remarks. Take $ for illustration and the rest 

will follow in likeness. 
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Then, it is quite clear that the residue contribution from the two parts 
are equal but opposite in sign to bring about a zero net contribution. 

Integrals dealt with up to the present have not involved elaborate 
treatment and most have been reduced to analytic simplicity, with the 
remainder purposely left incomplete for the moment. Next, attention will 
be placed on $ ^^^2 *^B12 expressed in Eqs. IV-16 & 17. 

First, recall the exponential integral is likewise expressible in terms of 
the incomplete Gamma function, i.e. E^(z) = - r (0, -z). r (ot,z) is an 
entire function of a, but in general, except when a is an integer, it is 
a many-valued function of z with a branch point at z = 0^^^. Since a = 0 
here, E^(-iK ^ F) is analytic in the ^ -plane. However, at ^ q, which is 
the zero of F, the argument vanishes and this causes the exponential 
integral to diverge. Thus, is a singular point. Expressed in series, 

« (-iK^F)" 

E. (-iK^ F) = Y+ In (iK^F) + (IV-29) 

where Y is Euler's constant. Clearly, in the limit as so that 
F 0, E^ presents there a logarithmic singularity which is defined as a 
branch point. For each K q when it is complex, there is also its complex 
conjugate Cq*. But since this point will always be enclosed by Cp^, only 
Cq needs to be concerned with. Now, aslEl"^ , iK ^F -> 1/^ 
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Certainly, E^(0) diverges, but note -iK^ FE. {-iK^ F) which for large 
behaves as 



By L'Hospital's rule 


1 im 

5-fco 


1 im 

^-KO 




00 
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Furthermore, 


lim ^iK^F ^ 


lim . 1 


and in ^nd •'^specti vely, 

g(g-M) . JL ^(C-M)L- _1 

Aj F ^ ^2 ’ A 2 F C 


Therefore, the integrands of ^/\22 *^B12 shown to behave as 1/^ 

or better for ^ becoming very large. This makes it possible to deform 
the integration contour in a manner identical to that illustrated by Fig. 
IV-8 used in the evaluation of branch point, contibu- 

tions arise only from those poles exterior of ^F+* as in the case before. 
Here, with the branch point at Eq, additional contribution will result 
from integration over the branch cut that must be made to account for the 
multi-valuedness. Thus, 
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e’ V E. (-iK^F) 



\ E.(-iK^F) dc 


(IV-30) 


where Cj^ is the contour around the branch cut which is to be discussed 
next. The character of Cq was examined in detail and portrayed by Fig. 
IV-10. Essentially, before time equals which is the minimum time for 
the waves to arrive at a given field point, stays on the real axis. 
After -^p, Cq moves onto ^FO to become complex. There is no unique way 
of making the branch cut but definitely it should be made in the most 
uncomplicating manner. For instance, it is highly desirable that the cut 
be straight and bypass t^FO and So, the cut can be made in the 
subsequent forms. Before t= -^p, depending on which side the point Cq 
may be positioned, it is made along the real axis away from CpQ as depicted 
in Fig. IV-11. This simplifies the integration in that it is performed 
with respect to u alone. Substitute the series expression Eq. IV-29 for 
the exponential integral into Eq. IV-30 and clearly, only the In term will 
contribute as the other two terms are single valued on Cj^. 


V 



Figure IV-11 Branch Cut Before i^p 


With the branch point not contributing over the vanishing circle 
t becomes 


'A12b 




K in (1. F) d5 

'^b/'^bb ^ ^ 


Let 


c _ iC(^-M) 
F “ “X 


K^e 


iK^F 


which is real along and But, 

In (iK^F) = In (c e^®) = In C + ie 


(IV-31) 


where C = |iK^ F| is the same but 6 is discontinuous on the two 
paths. For 0, 6 = tt on l^bo and e = - IT on and for > 0, 

6 = 0 on and 6 = 2 tt on This quite simply leads to: 
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Vl2b “ ^ 2 


Aii 

” / %<“> du 
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= i 2 n I Sp(u) du 

J^co 






(IV-32) 


Since the above is a purely imaginary quantity, it will not add on to the 
overall transmitted field. However, later in dealing with g, the same 
integral would yield a real quantity which will be shown to have a signifi- 
cance relating to the causality condition. 

The expression for ^/^i2b above is for t less than when 

^ is real. Now, turn attention to t greater than igp. Once more, the 
branch cut can be made straight and away from CpQ and Cp^ as depicted in 
Fig. IV-12 for the two possible cases of (1) Cq being complex and (2) 
^ Q being real but in between two Cpg's. 



Figure IV-12 Branch Cut After Tgp 
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The simplification of follows by employing again the 

procedure taken for t < tqp 
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where for t < t 
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iv) dv 


(IV-34a) 


with 


T ( c ) = K, e’’'?'' 


(IV-34b) 
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(IV-35) 


Note however that when j 3, -iK^F ®° . Express the exponential 

integral by its asymptotic expansion for large arguments, 


E^(-iK^F) 


e-iK^F 

iK^F 



and substitution into Eq. IV-35 leads to, for j = 3, 



The above shows that those points are no longer poles and hence reduce to 
zero. Thus Eq. IV-35 is replaced by 

^ B12s " e’ V E.(-iK^F) dE (IV-35a) 

C3 2 

This completes the analysis of "J* p with the exception of $^25 ^ g26 

which will be dealt with in the section on instability as they can be 
directly linked together. For the time being, an intermediate summary 
shall be made by collecting all the terms which have just been evaluated. 
The reason for having kept some terms in their incompleted forms becomes 
immediately clear as a multitude of cancellation takes place. So, write 


£1 



(IV-36) 
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Interpretation of these results will be reserved for a subsequent section 
but a brief description is in place here for clarity of the above equation. 
The first term (Eqs. IV-4 & 19) requires a simple analytic evaluation which 
is only complicated slightly in the determination of . The second and 
third terms defined by Eqs. IV-28 & 28a contain details on the resonance 
waves as indicated by the delta function operating on F{ ^ which is 
fairly simple to calculate. The definition of 4’/^i2b ^ B12b 

Eq. IV-31 to Eq. IV-34b depends on the time t and their evaluation require 
numerical integrations that may lead to some technical difficulties. 
Finally, the last two terms whose definitions can be found in Eqs. IV-11 & 
12 will be tackled later as mentioned. 


To complete the analysis, Eq. IV-16 for •I’g will be evaluated next. 
Since the course taken is identical to that for 4>p, the intermediary steps 
leading to the final solution will not be described. However, while there 
are similarities between ^ p and ^ q to high extent, there are also subtle 
but distinct differences and these are pointed out in careful examination. 
First, write down the results as follows: 

Qr.(0 

I 

(IV-37) 


4. =J- 

G ._3 




Re 


r 2 .^ 

d^ 




+ (J) 


D26 


^E26 


62 



where the expression for G is found in Eq. IV-4 and 
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(IV-37b) 


(IV-37C) 




W = K3e-¥<«3> . Qe{53) = ^ 

with defined way back in Eq. III-16a and in Eq. IV-28b 
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(IV-37f) 


It is immediately recognizable that there is a certain oddity about the 

first item in that this quantity can be non-zero before time equals t 

OG 

which is the minimum time for signals from the lower half of the ring 
source to arrive at a given field position. This appears to be in viola- 
tion of the causality condition established in Eq. IV-3, but a closer 
examination points to the resolution of this contradiction at once. It 
might be kept in mind that tfie causality condition simply requires 4> « to 
be zero before T-gg and not each term which makes up its expression. So, if 
the terms can be shown to sum up to zero, then the causality condition is 
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quite obviously satisfied. To do this, first examine the terms that make 
up <I>p where this problem is not encountered. The reason is there, in Eq. 
IV-36, the first term vanishes before t equals Tgp and furthermore, the 
terms O contribution as they are purely 

imaginary quantity. Now, returning to <I> g, it can be seen that while the 
first term is real before Tqq, $ ^^25 ^ E12b correspond to 

^ A12b il>g22b The fact that 

had been shown to satisfy the causality condition at the outset and then 
having properly carried out the integration all the way provides for a 
sufficient ground to assume that these terms have the self-cancellation 
quality. However, this point cannot be easily substantiated as direct 
proof is not possible. It will be necessary to perform the numerical 


integration on 5>Qi2b ^ E12b compare with the quantity of 
the first term. Next, observe that $ ^ 21^2 ^ E2K3 purely 
imaginary in contrast with the parallel terms in $p which are real and as 
mentioned can be related to resonance waves. Here, they simply do not 


contribute probably due to the influence of curvature effects. Lastly, it 


will just be remarked that $ ^25 ^ E26 again referred to as 
instability integrals and together with their counterparts in 4 > p will be 
treated in the next section. 
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V. THE INSTABILITY ANALYSIS 


This chapter is devoted to the examination of those terms identified 
as instability integrals in the last section. First, restate them: 


A26 


B26 




4— e'l'f dsdK 
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FjFjAj 
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^D26 is similar to ^ ^25 except in the absence of i and the function F 
being replaced by 6 . Taking away i and replacing L_ and F by L^ and G 
respectively in $^26 transforms it to the expression for $ ^ 25 * 
these integrals have in common is the fact that they all diverge. It can 
and will in fact be shown that this occurs only in a limited region in 
space and time as expected. Recall that Cg is a small circle of vanishing 
radius about the moving pole of 1/(K+ A 2 / ^2 '•s a 

function of K and by continuously varying K, the movement of the pole will 
follow a certain path C. This path C is merely an approximation to the 
actual path which is traced from the original eigenequation Sj defined 
in Eq. 11-24. In fact, may not be the only path of Sj. However, any 

other possible paths of Sj will not be be considered as they play no role 

in the stability analysis, which is the real concern here. That any path 

other than can be disregarded has been shown by Morgan^^^^ and Munt^^^^. 

The most central issue here is, how well does C approximate C^? Getting 
the answer to this question is not a trivial matter as remains an 





unknown because of its complexity which was the reason for making approxi- 
mations in the first place. So the alternative is to obtain not the exact 
but a fairly accurate picture by taking more terms in the asymptotic 
expansions. But before going on with that idea, examine first the curve C 
which can be constructed rather simply. Fig. V-1 shows a number of curves 
for various Mach numbers. Observe that in the limit of K tending to 
infinity, every curve independent of the Mach number moves toward the 
corresponding point of 1 + 0. This comes as no surprise since K 

equals infinite is the short wavelength limit where the plane vortex sheet 
analysis is valid and the zeros of 1 + F F g give rise to instability 
waves in that analysis. In the other extreme of K tending to 0, observe 
the curve moves toward and overshoots M slightly for the case of M = 0.5. 
This path is surely acceptable since the original path of integration 
can be placed above the curve in accordance with its requirement. But the 
figures show that for all the other Mach numbers, their corresponding 
curves C seem to move out toward infinity. If this was true, then would 
have been crossed in violation of the causality condition. However, 
considering the present problem to be well set up, there is no reason to 
believe that such incompatibility can take place and the paths are highly 
suspicious of being far from representing the true paths C^. The fact that 
C appears to fail in describing for small K is not entirely unexpected 
since the asymptotic expansions for large K were carried only to the second 
term. Now, one more term will be appended as suggested beforehand. Thus, 
the expansions for Hankel functions are: 
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with their arguments in between - ir 
functions are: 
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where these expansions are valid in the region - tt to tt . The expan- 
sions for the region of 0 to 2 tt deviate from the above in the manner 
consistent with the differences between Eqs. III-2 of the two-term expan- 
sions. Then, in reference to Eqs. III-4 & 5, the following definitions are 
derived. 
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Following some algebraic manipulations, and 02 in Eq. 1 1 1-8 now become 

Ol = - 1 (1 r FjFg) + gkr F,F, * ’’j) '’’l ' rj) ’ ^h^'> 

0 1 2 


i3 


(r. + r,) (r, - r,)^ (3 - 5 


o ^•‘ 0 / V'li •‘oj 

128(kr^rjr2) ‘ ^ y ^ 


O 2 = - (Fj - F^) (1 - FjFj) r gt,,V r, Cl ' '' 2 > O'! + T^Hl r 3F Fj) 

0 12 
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Without getting into the mathematical details, it can simply be stated that 
the above will transform the expression for the moving pole in the denom- 
inators of Eqs. III-15 & 16 to the form 
2 ^2 ^3 

+ ■ r, ^ - K + 


^^2^1 


(r^r2)^A^ 

where and defined as before in Eq. III-13a and 13b, and 


Ao = - ■ .L (fj + T2) (r^ - (3 - Er^r,) 


128 


1^2' 


Obviously, the quadratic equation will allow for two branch of curves 
traced by its zero as a function of K. In order to select the appropriate 
branch, first rewrite the above for its zero as 
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I6rjr2 U + T^T^) [ ^ 


6(3-5r^r2)(Hy2) pi 

(1 - 3 J j= 0 


(V-1) 


Then, it is at once clear that the plus branch is to be taken since near the 
zeros of expression reduces to 


K + 


i(Fj - r2)(l - 3Fjr2) 

8^2 (1 + rjF2) 


= 0 


which is dealt with when two-term expansions are used. Now, a new set of 
curves C are obtained numerically and sketched in Fig. V-2 based on the 
positive branch of Eq. V-1 and the same curves in Fig. V-1 have been 
superposed for the sake of comparison. For M = 0.5, the two curves are 
very nearly the same and approach the point ^ = M = 0.5 overshooting in 
the former and undershooting presently. For higher Mach numbers, the 
two curves are , decidedly different except near the K = “ limit where 
falling back to the points of the zeros of 1+ ^^^2 expected. In the 
case of M = 1.5, the curve appears to head toward the point C = M but then 
undershoots it drastically. For both M = 2.5 and 3.5, the curves again 
visibly diverge but not before the move towards their respective point of 
? = M has become more than evident. Once more, the reason for diverging 
in the direction of infinity is inescapably due to the smallness of K on 
that portion of the curves. It might be concluded then that the three-term 
expansion still falls short of yielding a representative picture of C^, 
thus requiring additional terms. On the other hand, there is no question 
about the notable improvement of the three-term curves over the two-term 
ones in snowing more sense of direction. Unmistakably visible is the fact 
that all curves now seem to move toward the point C = M, leading to the 
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speculation that is a smooth curve whose end points are at the zeros of 
1 + ^2 large K limit and at ? = M in the small K limit. 
Therefore, it would suffice to prove that the curve can indeed be traced to 
C = M in the limit of K ^0. This approach is preferred over manipu- 
lating on extra terms in the asymptotic expansions which can be rather 
cumbersome. For small K, it is necessary to exapnd Sj in ascending series 
and here to the second terms as follows; 



where y is the Euler's constant. Substituting these into the expression 
for Sj in Eq. 11-24 and collecting only the leading terms deliver the 
following equation for the zeros. 

j, (K r.)^ K r, 

Fj + T2C — In (~Y^) = 0 (V-2) 


Clearly, when K = 0 the above reduces to a zero at C = M and on first 
inspection it appears to be a zero of order two. But noting ( E - M) is as 
well a part of the numerator in Eq. 11-23 makes E = M just a simple pole. 
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thus in consistency with the poles for all values of K. This completes 
the proof on where begins its path at the K = 0 end. Now Eq. V-2 can be 
used to show, as illustrated in Fig. V-3, how it approximates for small 
K. Note the curves from Figs. V-1 & 2 have been superposed so that the 
combination of all three sets of curves can help to shed light on the 
general nature of C^. Here, those curves resulting from two-term expan- 
sions are labeled Cp Cg for three-term, and for small K. For M = 0.5, 
Cj is sandwiched in between and C 2 in the most conspicuous manner that 

Cj can be safely assumed to lie within that narrow strip, beginning at 
M and ending at E = E p the corresponding zero of 1 + Fg which is 
displayed in Eq. IV-2. For M = 1.5, there is no doubt in seeing to 
approach the other end point at E = C p so a dashed line has been 
extended from to complete a curve most likely to portray C^. For either 
M = 2.5 and 3 . 5 , seems to head a little off course but still appears to 
proceed toward the other end, so a dashed line is drawn in each case to 
simulate the likely profile of C^- 

It is established that the path C traced by the two-term expansions does 
not come close to approximating C^, especially in the small K limit. The 
situation is remedied by taking one extra term and at the same time 
approaching from the small K end by using ascending series expansions. The 
results have been shown to furnish sufficient grounds for drawing up a 
conclusive picture of C^'s profile. As remarked in the beginning of this 
section, has not been found exactly due to the complex nature of Sj, but 
the approximate picture obtained is all that is necessary in showing the 
moving pole to be bounded. So in working with ^26 
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FIGURE V-3 M0VIN6 POLE PROFILE 
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FROM SERIES EXPANSIONS 








taken to move along in that fashion. With this knowledge, the instability 
waves can be examined for the spatial extent in which they propagate. The 
moving pole illustrates specifically the presence of a singularity for 


every value of K, where the one corresponding to K = ® is precisely the 
singularity found in the line source problem that gives rise to instability 
waves. The integrals under consideration, i.e. ^>^25 etc., can therefore 
be regarded as made up of infinitely many integrals, with each given rise 
to instability waves at a distinct K. The idea then is to get a picture of 
how the region of propagation varies with K. First, express the integrals 
in a representative form. 


1 

471 \ 5 - £5 


(V-3) 


with continuously varying for all K such that for the integral corre- 
sponding to K - CO , or ^ ^2 ^ 5 ®ce no longer included 
since ^ 2 always in the lower half while ^ g moves into the lower half 
as soon as K # eo and it was shown in the last section that any poles below 
the real axis can be enclosed within Cp^ to nullify their contributions. 
Inversion of the above can be shown to lead to Qj( ^g) • 6 [f(c 5 )] • 
Qj(^ 5 ) represents the magnitude which is generally expected to be large. 
For M< 2 it is supplemented by (see Eq. IV-5a) the Heaviside 
factor containing the information on when is exterior of Cp^ which is 
directly related to space and time. Observe the delta function has a 
complex argument. Using the notation of Eq. IV-5, that is for F( ^g) = “ 

+ i 6 it can be shown that 
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S(a+iB) = f - 6<">(a) 

Thus, the instability waves propagate along a = 0 in the region 3 < 0. 
The spatial extent can be determined by eliminating time between the two to 
find it boundary. Recall first 

F = Z + Tj (R - 1) + T2 (1 - R5) - Ti 

ForK =“, = M/2 + in and 

a = Z + I Vm^ + 1 -ij (R + R3 - 2) - 

B =n[^ C + 1 + 1) (R - Rs) - t] 

from which the boundary line for the instability waves becomes 

Z = 4“ [ + 1 + 1] (R -R5) - -1] (R+R5-2) (V-4) 

Since dZ/dR = 1, it is clear that the line makes a 45° angle with the axial 
direction. For K , the boundary equation becomes, with ^ g(K) = u + 

iv, 

z = ^ [fiC' - 1) ^ * r,(l . Rj)] 
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and 


dz 

dR 


I r, - Rer 
V mi 


This means the boundary no longer makes a 45° angle but still is a straight 
line. The expressions of Cg for K ~ are not known explicitly and must 
be resorted to approximations provided by Fig. V-3. By tracing from K = ® 
along the possible profiles for Cg, the values of Cg can be roughly 
estimated to be used in Eq. V-5. For illustrative purpose, the boundary 
lines thus obtained are shown in Fig. V-4a for the case M = 0.5 since Cg is 
believed to be most accurately predicted there. The boundary lines for the 
instability waves from both 4>p coming first and 4> g second are included. 
It is clear from these figures that the boundary line for K = ® is always 
the outermost while the other lines, for K moving away from » , all fall 
to the right with lesser angles. The instability waves propagate in these 
regions where B (K) < 0 as singular wave forms defined by a(K) = 0, 
which are straight lines more or less perpendicular to the boundary lines 
as shown. At this point, it might be recalled that as K becomes smaller, 
the corresponding point on Cg approaches the point ? = M. What this 
means is that because C = M is always enclosed by ^FO’ somewhere along the 
way as K becomes small, say at some critical value K = Kg, Eg will cease to 
be outside of CpQ, which is equivalent to being enclosed by Cp^ as shown in 
the last section, and hence for K< Kg, instability waves may no longer 
appear. This relationship is illustrated by Fig. V-5. It must be 
recognized that the present analysis is only an approximation, and just as 
the explicit expression for E g as a function of K remains undetermined. 
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V 


Decreasing K 



Figure V-5 The moving instability pole and its connection to the 

region of propagation 


the value of and consequently that of are approximation, and just as 
the explicit expression for C g as a function of K remains undetermined, 
the value of and consequently that of C q are unknown quantities. 
Therefore, the precise location of the boundary line for K = cannot be 
established. But quantitatively, is of the order of 1 and since K = kR^, 
Kq is of the order of the jet radius. The possibility remains for that 
line to come in contact with the vortex layer. Then the instability waves 
may be characterized as having their presence in a limited region which is 
a function of wavelengths such that for short waves, the boundary makes a 
45° angle with the vortex layer and as wavelengths become large, the 
boundary moves to the right and draws near the discontinuous surface until 
finally, at some critical wavelengths, the boundary comes to fall on the 
vortex layer to end the instability propagation for any longer waves. This 
characterization is merely an assumption and it is actually more 
defensible to simply state that the instability waves have wavelengths on 


80 



the shorter end. Even though the profile of for M = 1.5 is less 
perceptible from Fig. V-3, an attempt has been made at obtaining the 
corresponding boundary lines and these are shown in Fig. V-4b. In contrast 
to the lines for M = 0.5, these are tighter together and the angular 
variation is nearly absent. Similar pictures are expected for higher Mach 
numbers so the previous contention that the boundary line for K = may 
fall on the vortex layer is less probable. 


The preceding discussion has centered on the case of M < 2 2. 
Basically, the curvature effects are responsible for the predicted func- 
tional dependence of the spatial extent of instability waves on the wave- 
lengths. Next, consider what happens for M > 2 2. In the plane vortex 
sheet model, the instability waves will cease to exist as the complex poles 
become real to characterize the neutral stability waves. However, an 
examination of Fig. V-3, for the case of M = 3.5, reveals the vortex layer 
model to yield a different picture. The profile of Cg clearly shows that 
the pole is truly on the real axis only when K = ® and as K is varied from 
large to small, it moves into the complex domain. Since in reality, K = ® 
corresponding to zero wavelength is just a limit, the overall picture calls 
for the continued presence of instability waves at the higher Mach numbers. 
The representative Eq. V-3 can be used once again to ascertain the spatial 
extent. Here, however, when K = " , the inverse is accompanied by H^ 
defined in Eq. IV-7. Eliminating time between this expression and F( = 
0, which is now real, yields the boundary line equation as follows. 


Z 





(I - Rj) 
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where and are given in Eq. IV-2. dZ/dR is less than unity so that 
the angle which this line makes with the horizontal is less than 45°. When 
K ^ ® becomes complex and the boundary lines can be determined 
identical to the M< 2V^case. These lines as expressed by Eq. V-5 are 
again expected to fall to the right and closer to the vortex layer until at 
some critical wavelength, the pole moves inside of CpQ to end the insta- 
bility wave presence. 

In summary, it has been shown the instability integrals will give rise 
to instability waves for all Mach numbers but that they are wavelength 
dependent and cease to exist for wavelengths large compared to the jet 
radius. Furthermore, the region of propagation is similarly affected in 
that it becomes less extensive for longer wavelengths. 
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VI. RESULTS AND DISCUSSIONS 


To validate and interpret the results obtained in the previous 
chapters, a computer program based on the analytical solution was 
developed and numerous cases have been investigated. These numerical 
results are now presented and discussed in this chapter. 


First, the organization and method chosen for analysis are summa- 
rized. In Chapter II, the problem is formulated and its solution in the 
Fourier transformed space is obtained in the form of a multiple integral 
whose integrand involves complicated Bessel and Hankel functions. Chapter 
III contains the derivation of an approximate version of the integral 
solution for the transmitted field; this is based on the assumption that 
the jet and source radii are sufficiently large so that the two leading 
terms from the appropriate asymptotic (large argument) series expansions 
of Bessel and Hankel functions can be used. In approximating with these 
expansion terms, contributions arising from internal reflections are 
isolated and neglected. Singularities of the integrand are examined in 
Chapter IV prior to the inverse analysis which is carried out in the spirit 
of Cagniard's technique of contour deformation. 


The results are now 


<I> = ^)p + 


477 

1 

477 ^ 1 ^ 5 ^ 


Re 

Re 


'SF 




put together and presented below. 
^A2K3 * "^BBRS * ^A12 ^B12 ^A26 

^D12 ^E12 ^D26 ^EBbI 
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where ^ p and $ q represent the fields from upper and lower half of the 
source respectively. Each is subdivided into component parts as 
specularly transmitted, instability, and neutral stability waves. The 
instability waves ^ B26 ^ E26^ requiring a separate analysis 
are discussed in detail in Chapter IV. And while both halves of the source 
may give rise to instability, the neutral stability waves are caused only 
by the upper half. They are, like the instability waves, confined to a 


limited wedge-like region whose boundary is described by the elimination 

, r X M(R - R.) 

of time between F(5^) =0 and t - , ^ — =0. A field that would 

VM^-4 

modify the specularly transmitted field in the region of its propagation 


appears to follow these neutral stability waves, as represented by 'I’ ;^2K3 
and $g2K3 IV-28. 


The remaining terms in Eq. VI-1 form the specularly transmitted 
field. If the solution were approximated by just the one-term asymptotic 
expansion of Bessel functions, ‘I' would have been absent because its 
prediction is a consequence of accounting for the finite curvature which is 
lacking in the simpler approximation. Since upper and lower halves of the 
ring source have their curvatures opposite in sense, it is not surprising 
to find 4>gp and somewhat different in characteristics as pointed out 
previously. In fact, every term in $ p is different from the corresponding 
one in 4’g, although the mathematical differences would not be meaningful 
unless they can be interpreted and explained from the physical point of 
view. 4>gp and given in Eqs. IV-19 & 37 are fundamentally unequal in 
the order of arrival at a given field point. Arising from the upper half 
closer to the point, $gp clearly arrives ahead of <t> gg. This means that 
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<I>jP should have the identity of a wavefront while 4>jg, which comes in 
succession, should not because the total field 4> cannot possess but a 
single wavefront. As mentioned formerly, the only way to substantiate the 
above is by working out numerical evaluations and some examples will be 
presented a little later. 


Next, attention will be focused on the last remaining terms of $ 

4>bi2* ^012 ^£12 IV-32 to 34b and Eqs. IV-27d and 

37e. These integrals cannot be simplified analytically any further as 
shown in Section IV and again must resort to numerical evaluations. A 
significant point to look for is the way in which they will modify and 
to form the total specularly transmitted field. There is no 

uncertainty in p to satisfy the causality condition as both 4 )gp and 
and $ s**® ^ero before the time of arrival. But the same is not obvious 

for 4>g, the field coming continuously after the wavefront. It is expected 
that $ gg and <3 >qj^ 2 ^ £12 ^ value for the satisfaction of its 

causality since the two parts are non-zero separately. The numerical 

evaluations, which will be discussed now, indeed show this to be the case. 


The calculations of $ gp and ^ gg are fairly straight forward because 
of their analytic simplicity. For a given set of field and time data, it 
is just a matter of finding the zeros of F and G for substitution into the 
expressions. The minimum time of travel or the time of expected wavefront 
arrival for (j> is manifest in $p. This time is directly related to F in 
that it corresponds to the time when Re(F) first intersects Im(F) at one 
point of CpQ on the real axis. One can find this time by trial and error 
numerically but to do so is both costly and unnecessary. For the 
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specularly transmitted field is something that is theoretically predict- 
able from the approach of geometric acoustics. Surely the complex nature 
of the present problem is sufficient to discourage any such attempt but at 
least one information can be readily extracted without much difficulty. 
And that is precisely the minimum travel time. First, the minimum time for 
the signal to reach the surface of the jet is given by 

^ ^ ■ h 

In this time, the source is imaged downstream to the point 


And the minimum time to reach the field point (Z,R) strictly evaluated by 
geometric acoustics is then given by 

For Z < Z|p 

~0F " ’'^IF + (R - 1)^ M >1 (VI-2a) 

= T2p + + (R - 1)^ M <1 (VI-2b) 


where 


= z - Mt2p ± ( 


2F 


' IF) 


86 



and + is for Z > Zj^p. For Z > Z p, T^p is the same as Eq. VI-2b for 
both M less and greater than 1. In the above, ^p and should yield a 

minimum in XQp. Therefore, digp /df 2 p = 0 or 

^[Et " V ^ - E^ [m (t|p - T^p)^ 5 T2p)]= 0 

Tqp is the minimum time of travel for $ p as well as $ so that a wavefront 
is expected then. On the other hand, tqq is just the minimum travel time 
for <I>g and it is different from xgp only in the expression for which 
is now (I+R 5 ) instead of (1-Rg). In the examples to follow, x Qp and x qq 
are calculated by the above expressions. 


Numerical integrations of difficult to perform. 
Considering the semi-infinite interval of integration, an initial attempt 
by Laguerre fntergration proved unsuccessful because the integrands do not 
decay in a simple exponential fashion. More often than not, the path of 
integration would pass relatively close to a singularity causing the 
integrands to fluctuate. Therefore, the task is completed by Gaussian 
integration instead at the cost of considerably increased computation 
time. This involves subdividing the semi-infinite span into small inter- 
vals and integrating with Gaussian weights from the starting point. Each 
subinterval is integrated until the desired accuracy is attained. Then, 
the same is performed onward until two intervals in sequence do not differ 
in value greater than a specified amount. If the maximum number of 
Gaussian points used do not lead to a satisfactory result for any sub- 
interval, then it is subdivided even further repetitively until necessary. 
Since the integrands decay smoothly and rapidly for large arguments, that 
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is, when the path leaves the region of singularities, convergence to three 
place accuracy usually occurs well before the argument reaches the order of 
50000. Convergence is slower when ris large because the integrated values 
are much reduced. The results presented will be divided into two regions. 
They are designated as the quiet and the noisy regions depending on the 
absence or the presence of instability waves. In the quiet region, inte- 
gration always converges. But in the noisy region, convergence is 
difficult to obtain once entering the actual regions where instability 
waves propagate. This is because the branch cut integration now passes 
through the moving pole contour which influences the integrand to 
oscillate with large amplitudes. In reality, even if convergence is 
achieved, the transmited field thus calculated still does not include the 
contribution arising from instability waves. In order to predict the 
amplitudes for the instability waves, non-linear effects such as due to 
viscosity must be retained in the conservation equations. The present 
analysis based on the linear theory can only give a qualitative picture of 
where and when instability waves propagate and how large a region in space 
and time they occupy. Thus in the results to be presented, the instability 
regions will be indicated simply by blanks in that portion of the curves. 

From the solution for the velocity potential, two field variables of 
more interest can be derived, namely the acoustic pressure and the vortex 
layer displacement. In the linearlized Euler equation, pressure is 
related to velocity potential by its time derivative which can be evaluated 
in analytic simplicity as shown in Appendix B. The derivation of vortex 
layer displacement found in Appendix C is more involved in that it requires 
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the same inverse analysis applied to velocity potential. Numerical calcu- 
lations for all three field variables are performed by the computer program 
listed in Appendix D. 

First, analyze what the solution predicts in the quiet region. Figs. 
VI-1 through 6 are plots of velocity potential $ vs time for various 
field points, Mach numbers, and source radii. There are three curves in 
each plot. The dotted ones are p, the dashed ones and the solid 
ones ^ which is the sum of the other two. Several general features are 
immediately evident from these figures, some of which are rather signifi- 
cant. First of all, 4>p admits a wavefront at the time of its arrival as 
expected. This is represented by an acute rise from zero at that instant 
and then followed by a sharp fall before a smooth decay takes place. Next, 
observe the behavior of ‘J’g at its arrival which comes slightly behind $p. 
Numerically, the sum of and ^ ^£12 shown to add up to zero 
within the desired accuracy. Thus, the causality condition is satisfied in 
the fact that3>Q rises from zero. The total specularly transmitted field 
hence possesses just a single wavefront. These figures demonstrate the 
fundamental differences between $p and To describe the sound field 
for the same data points, now look at the pressure vs time plots, which are 
presented in Figs. VI-7 through 12. The first five of these are for small 
ring radius, = 0.1, so that the source is nearly centered on the jet 
axis. Fig. VI-7 is for a point downstream of the source fairly close to 
the jet. There isn't much difference except in the amplitude between it 
and Fig. VI-8 which is at the same downstream point but farther away from 
the jet. The reduced amplitudes and the smoother decays mean lower sound 
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pressure. If now the point is moved further downstream while keeping the 
same radial distance as in Fig. VI-9, a significant change in the curves 
are clearly observed for supersonic cases. Here the amplitudes are much 
increased and exceed the corresponding cases in Figs. VI-7 & 8 which are 
both for points closer to the jet. It is clear from this that more sound is 
refracted in the downstream direction. What this means is that with 
increasing jet speed, more sound is refracted downstream. This is even 
more evident from Fig. VI-11 which is for a point straight out from the 
jet, or more upstream than the other cases. While its distance from the 
jet is no greater than the point in Fig. VI-7, the amplitudes are much 
reduced indicating that very little sound is transmitted upstream, 
especially in the supersonic cases. Between Figs. VI-8 & 9, one can detect 
the effect, but otherwise the amplitudes are pretty much dependent on 
distance between the jet and the field point. The cases discussed thus far 
are more or less representative of the near field characteristics. If the 
point in Fig. VI-9 is moved outwards in the radial direction, but fixed 
axially as in Fig. VI-10, then it is out in the far field. A comparison 
reveals some interesting features. In addition to the amplitude loss due 
to radial decay, there is the apparent stretching in the time decay history 
so that it takes longer period for the quieter sound to fade away. The far 
field radiation also seems to affect the transmission characteristics for 
supersonic jets. The curve for M = 1.5 in this case is no longer distin- 
guishable from that for the subsonic speed, and for M = 2.5, the curve 
appears to start losing its supersonic identity. Finally, the effect of 
source radius can somewhat be examined in Fig. VI-12 which is for the same 
field point as in Fig. VI-7. As the radius increases, the transmitted 
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field takes on a different look from previous cases. The relative lower 
amplitudes reflect the influence of the square root quantity on source 
radius in the denominator of the solution. The character of the pressure 
disturbance is also affected due to the fact that as source radius becomes 
larger, more time separates the signals from various parts of the source. 

Next, examine the disturbed fields in the noisy region closer to the 
jet axis. In this region, vortex layer displacements as well as velocity 
potential and pressure can be examined. As mentioned earlier, instability 
regions are the blanked portion of the perturbation curves. But as the 
exact time of passage in the K = ® limit can be calculated, that time is 
indicated by a short vertical line in each case. In Fig. VI-13, the field 
point is relative downstream making a 30° angle with the jet. Velocity 
potential is plotted on the top and pressure on the bottom for two Mach 
numbers. As expected, more sound is refracted into this region. While the 
disturbances decay rather smoothly in the quiet region, here they actually 
become increasingly unstable until the passage of instability waves is 
observed. After its passage, however, the disturbances die out quickly so 
that the noise covers just a short time period. This is especially true in 
the supersonic case where the impulsive nature of the source is very 
evident. Fig. VI-14 is a plot for vortex layer displacement vs time at two 
field points. On the upper left corner, the point is directly in line with 
the source relative upstream. For subsonic flows, represented by M = 0.5, 
the displacement is small and smoothly decays. For supersonic flows, as 
disturbances cannot propagate upstream, the vortex layer remains undis- 
placed for all times. The other three curves are for a point somewhat 
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downstream. In all cases, the vortex layer becomes excessively displaced 
and finally ends up with instabilities. But after its passage, vortex 
layer returns to quiet very soon. For M = 2.5, besides the instability 
disturbances, there is also the resonance disturbance which gives rise to 
the resonance modified field. This is indicated by the increased displace- 
ment following the passage of instability waves. This field which follows 
the resonance wave is shown to die out quickly. It may be concluded that 
the solution not only predicts a region of instability but also a region of 
resonance. Finally, Fig. VI-15 shows the time lapse downstream direction. 
It can be seen that for small time a disturbance begins to develop. This 
includes a wavefront followed by the instability and a smooth decay. Then 
as time increases, the disturbances move downstream with increasing ampli- 
tudes behind the wavefront. This large displacement is always quieted very 
quickly after the passage of instability disturbances, which is in agree- 
ment with the previous two figures. 

In conclusion, an approximate solution to the problem of transmission 
of sound across a cylindrical vortex layer has been obtained. The results 
are considerably different from the plane vortex sheet case because of the 
added role played by the curvature of the jet. In comparison with the 
plane case, the specularly transmitted waves are more complex and requires 
some numerical integration. Resonance waves are identically predicted for 
M > 2, but there is also a wave field whose modified effect appears to 
extend the region of resonance just as the instability waves cover a region 
in space and time. The instability waves are predicted to exist for all 
Mach numbers but vanish for wavelengths large compared to the jet radius. 
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And the region of propagation is similarly wavelength dependent. The 
prediction that large wavelengths are unimportant in the transmitted field 
is in agreement with some experimental observations. 
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Legend to Figures VI-1 to 15 
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Non-dimensional quantities are defined as follows: 
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APPENDIX A. SINGULARITIES OF THE INTEGRAND 


a. The zeros of Aj = ( + r2)(l + Fj F2) have been determined in the 
plane vortex sheet problems as is evident in Eq. III- 9 , but they will be 
verified here for the sake of completeness. Recall 

Now let Fj = Pj e ^*^1 
F2 = P2 e = a^ + ib^ 


where bj^ and b2 have the same sign. Then, Fj + F 2 = 0 requires a^ + 82 = 0 
and b^ + b2 = 0 . The only possibility for satisfaction of the latter is if 
“ ^2 ” translates to Pjsin92 = - P2sinP2, it requires 

either 6^ = 0 and 62 = n or vice versa. The condition a^ = -a2 means 
p j^cosSj^ = - P2COS0 2 and it is clear from Fig. IV -1 that 6^ = 0 

and 0 2 = 7 ^ while P^ " P 2 “ M > 2 . Thus, F ^ + F2 has a 

single real zero at C = ^3 = M /2 for M> 2 . 







+ 0 


22 


2 


jjn 



1 + = 0 requires a^a2 - b^b2 + 1 = 0 and a^b2 + a2b^ = 0, from which 

cos(e^ + 82) =- — and sin(6j^ + 62) = 0 are the consequences, leading 

to + 02 = 71 and ~ ^’’S* ^V-1, it is obvious that = tt 

only 1f Ur, the real part of the zero C„ . „r + ,VR, equals M/2. Then, 

Pi = P2 ■ ^ imaginary part of can be determined. 

Pi = Pii P12 = Vvr 2 f (-^- 1)^ 1 /vr^ + (-^+1)2 = 1 

Square both sides and solve for v. 

Vr^ = - (M^ + 1) ± /m^ + 1 


The spurious negative root is discarded. 

,2 


-[/ 


+ 1 


- (-^+lr 


+ 1 ) - 


Vm^ + 


H. 
A 


Then, 

for M < 2-^ 
for M > aVi" 


Thus, 1 + riF2 have two zeros, which for M< 2V^are complex conjugates and 
for M >7/2 are real. At M = 2 ^, Vr = 0 and the two zeros collapse into a 
single value of M/2 which coincides with the zero of r 1 +T2* 

b. The zeros of A2 = i(ri^ - F2^ ) (1 - 3rir2)/8. The first part of 
2 2 

"^2 ~ ^ ^2^ ■ ^2^ examined in a). By making use of Eq. IV-1 
again, ^1 2 “ ^ requires PiCOsQ 1 = P2COSO2 and PiSin6i =P2Sin^ and the 
only possibility of this is if 6i = 62 = tt /2 and Pi = P2 when M< 2. The 
corresponding zero is then at M/2 as before, which is a point on the 
overlapped portion of the branch cuts. 
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Any zero of 1 - must satisfy the conditions 1 + 3(b2b2 - 2222^ “ ^ 

and a|b2 = which translates to cos e^sin^ = -sin02cos62 or 

sinCe^ + 62) = 0 . Since b^ and b2 have the same sign as required by the way 
branch cuts are made, a^ and a2 similarly must have the same sign as well. 

This leads toe j ■*‘62 *^ 1*^2 " " which is not possible under any 

circumstances and therefore the conclusion is that 1 - 3F2T2 
any zeros. 

c. The zeros of and L_ are similar to the zeros of F^ +F ^ andF ^ ■ ^2 
but their presence and positions will depend on and R. As R tends to 
become much larger than R^, the zero of will move very closely to 0 on 
the branch cut of F^, while that of L_ is likely to disappear. 

d. The zeros that have been identifed thus far all fall in the C-plane 
and are functions of only the Mach number. The zeros of K + A 2^^^! ^2 

on the other hand are identifiable in both the f;- and K-planes in addition 
to being Mach number dependent. Any zero in either plane is a function of 
the other variable and the resulting movement makes it rather difficult to 
handle. In the K-plane it is straightforward to recognize the zero as 
function of or K = = - A2/(F2F2 A^). To find the corresponding 

zero as function of K in the C-plane is more elaborate and will require 
search by numerical means except for the limiting cases of K 0 and 
K . However, as these are directly linked with instability, they will 
be treated in more details later. 
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The zeros of F and G have different physical interpretations so their 
treatment is given in the main text. In summary, all the zeros of Eqs. 
III-15 & 16 have been identified or pointed out. Those which are poles to 
the integrand will be shown to contribute to instability and resonance 
associated with the vortex layer model. 
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APPENDIX B. DERIVATION OF PRESSURE 


From the linearlized momentum equation, the acoustic pressure p is 
related to the velocity potential by 

P = - P (B-l) 


In terms of non-dimensional quantities, this relationship reads. 


P = 


^0 P 

pa^Q 


9t 


(B-2) 


where and t are defined in Section III. 


Here, Eq*. A-2 will be calculated for the specularly transmitted wave 
field. As ^ ^G* total pressure field will also be the sum of 
Pp and Pg. This is expected to introduce some errors at the vicinity of 
arrival of Pg since the derivative of $ g is discontinuous there. Only 
the derivation for P p will be demonstrated with the understanding that 
similar procedure applies to Pg. First write. 



Now, Eq is such that, from Eq. III-16 
F (Cj,) = Z ♦ rj(£^) (R - 1) + rj(£^) (1 - Rj) - = 0 
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and so 


T = '’ 2 < 5 o> - "s>] 


Thus, 

3t _ ^ 


r C. (R-1) (5o-M)(l-R0 T r 1 

[ r° (g ]-[^^ rj(c„)(R-l) . r^U^Xi-Rj)] 


(B-4) 


From Eq. IV-19 


Qf (C„) 


5o<5o 


where 


Ai(^o^ = (Tj + T2) (1 + rjF2) 
Thus, 


C=C 


0 


A,(f )(2F -M) - F (r -M) A' (E ) 
— Qf(?„) = Qf(5„) = 5 i-5- 




(c,) 


(B-5) 


with 


2[V2'«o> * <«o-”>ri(5o>] " «0<V”>[i7(y 


(5o-”> . 5o 


w 


(B-5a) 


Next, 


dF 

di 





€q(R-1) (^o-M)(1 - R5) 


(B-6) 
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Then, 




Co (R-1) 

rj(Cp)(R-i) - lYTc^l 

^ 


(Co-M)^U-R.) 

■n 2/_ \ 3f. 


r/(s„) 


(B-6a) 


Eqs. A-5 through 6a apply to the first term 

, T F'(c„)Q^(c,)-Qp(gr(g 

SE L ^ j —rm 


’0 dC 


F'(Co) 


The second term of Eq. A-3 is made up of the sum of the time derivative of 


the two integrals ^’aiou and <J)| 


As the two are identical in nature, 


it suffices to show the derivation of just the first part. From Eq. IV- 


^A12b = " 


■/: ", 


(u^ + iv) dv 


where 


Sp (Uq + iv) = S^(C) = e 


and F is expressed in Eq. III-16. By Leibnitz's rule of differentiation. 


^^A12b 

3t 




3 Sp(uQ + iv) 


(B-7) 
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H \ Sp (Uq + iv) 


(B-7a) 


(B-7b) 
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APPENDIX C. DERIVATION OF VORTEX LAYER DISPLACEMENT 

The vortex layer displacement is related to the velocity potential 
through the kinematic boundary condition applied at the interface r = Tq, 



dr at 

Expressed in the Fourier transformed space, 




118 



where 


rv 1 / 1 ^ 3 V . 

8k ^ r I’j, ' r„ r ^ 

St 0 1 

F - 1 -L < ^ 3 ' 

1 - 1 ■ 8k ' r^Fj > 



The above equation is then evaluated in the exact manner as described in 
Sections IV and VI. 


119 



APPENDIX D 


C « K K « M » M K II KM « mOI M )( N K NNX K K K « M K IM» M KK 

C MAIN PROGRAM 
C 

C R - RADIAL DISTANCE 

C RM - MACH NUMBER 

C B - RING SOURCE RADIUS. 

C A - TIME 

C X “ DISTANCE IN THE DOWNSTREAM DIRECTION 

C ISP - 1 BEFORE WAVE ARRIVAL. INITIAL INPUT 

C - 2 AFTER WAVE ARRIVAL. OUTPUT FROM ROOT FOR COMPLEX ROOT 

C - 3 SAME AS 2 BUT FOR REAL ROOT 

C NI • MAXIMUM NUMBER OF INTEGRATION SUBINTERVALS. 

C 01 - ARRAY FOR ESTABLISHING SUBINTERVALS. 

C NL - NUMBER OF FUNCTION EVALUATIONS OF FUNCA OR B. 

C NG - NUMBER OF SUDINTERVALS NECESSARY TO ARRIVE AT 

C THE FINAL VALUE. 

C NT - NUMBER OF TIME INPUTS 

C MAXIMUM OF 30. STOPS EXECUTION IF 0 

C NC - CASE NUMBER 

C NW - 1 FOR VELOCITY POTENTIAL 

C - 2 FOR PRESSURE 

C - 3 FOR VORTEX LAYER DISPLACEMENT 

C FCR « 0 FOR PHIA, -1 FOR PHIB 

C FCI - -t FOR PHIA, 0 FOR PHIB 

C TIME - SUBROUTINE FOR GENERATING TIME INPUTS 

C ROOT - SUBROUTINE FOR FINDING THE ZERO OF F OR G 

C SHOD - SUBROUTINE FOR CALC’ULATI^IG MODIFIED FIELD IN THE 

C NEUTRAL INSTABILITY REGION, WHEN M > 2. FM ITS VALUE 

C SPEC - SUBROUTIt^E FOR CALCULATING THE SPECULARLY 

C TRANSMITTED FIELD, AND FA ITS VALUE 

C LAG - SUBROUTINE WHICH PERFORMS GAUSS 

C INTEGRATIONS ON FUNCA OR FUNCB, AND FCB ITS VALUE. 

C F - THE TOTAL FIELD, PHI. 

C 

IMPLICIT REAL**6(A-H,L,0-Z) 

DIMENSION T(30),DI(50) 

COMMON /CMA/ RM,R,B,A,X /CM3/ U5P,FCR,FCI /CMC/ NW 
EXTERNAL FUNCA, FUNCB 
EPS=1 .ODO 
5 EP5=EPS/2.0D0 
T0L1=1 .ODO^EPS 
1F(T0L1 .GT.1 .ODO)GO TO 5 
PI=^#.ODO<^DATAN( 1 .000 ) 

READ 10,NI,NW 
10 FORMAT(212) 

READ 15, (OKI), 1 = 1, Nil 
15 FORMAT(7D10.2I 
20 READ 25,NC,NT,X,R,B,RM 
25 FORMAT(2I2,^D10.2) 

IF(NT.EQ.O)5TOP 

NL = 0 

NG=0 

RMH=RM»»0.5D0 

DF=2.0D0^PI*»D5QRT(D^P) 
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CALL TI^ECT•EP5,^^^) 

PRINT 30,NC,RM,B,R,X 

30 FORMAK '-S'CASE M= MPD12.5,*. B= *^D\Z.Sp\ R= **012.5. 

♦ •* X= •.012.5//5X,*TiriE',nx**U**11X,’V*10X,*ADS’.CX**SPECSeX. 

♦ •NODI* *12X, •INTO* .nx, ‘PHIA* *8X, *PHIB* *9X, ’PHI V) 

I5A=1 

ISP=1 

DO 65 1=1, NT 
FCR=O.ODO 
FCI=-1.CD0 
A=T(I) 

1SB=ISP 

1SP=ISA 

DO 70 J=l,2 

CALL ROOT(U,V,ISP,FS) 

USP=U 

FCB=O*0D0 

FA=0.0D0 

Fn=0.CD0 

1F(B.LT.O.ODO)GO TO ^i5 
Ft1=SnOD(RMH) 

GO TO (60, AS, 55), ISP 
A5 FA=5PEC(U,V)/PI 

GO TO (50, 55, 55), ISP 

50 FCB=LAG(FUNCA,NL,NG,U,HI,DI,ISF»)/PI/e.ODO 
GO TO 60 

55 FCD=LAG(FUN’CB,NL,NG,V, HI. Dl, ISP )/Pl/C. ODD 
60 F=(FCD+rA+FM)/DF 
B=-B 

CO TO (65,80),J 
65 FT=F 

FCR=-1 .ODO 
FCI=0.0D0 
ISA=I5P 
ISP=IS3 

70 PRINT 75,A,U,V,FS,FA,Fri,NL,FCB,NG,F 

75 FORHAK • • , 1 PDl 2 .5 ,5( 2X,D1 0 .3 ) ,2X.I3, 1X,01 0 .3, 1X,I2,2X,Ot 0.3) 

80 FT=F4FT 

85 PRINT 90,U,V,FS.FA,HL,FCD,NG,F,FT 
90 FOPMATC •,1AX,A(1PD10.3.2X),12X, 
♦I3.1X,D10.3,1X,I2.t2X,2(2X,D10.3) ) 

GO TO 20 
tUD 

SUBROUTINE TINE( T.EPS.NT ) 

C 

C PROGRAM FIUOS: MHEU HT=30, THE MIIIIMUM TIME OF TRAVEL FROM 

C SOURCE TO THE FIELD POINT, AND GENERATES 30 TItlE STEPS FOR INPUT. 

C WHEN NT < 30, READS IN NT TIME STEPS FOR INPUT. 

C ZERO - SUDROUTINE WHICH SOLVES FOR THE MINIMUM TIME 

C AMM - MINIMUM TIME 
C 

IMPLICIT REAL?‘3(A-ll,0-Z) 

DIMENSION T(30) 
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COMMON /CMA/ Rtt,R,B,A,X /CTA/ AZ*S,R2 

EXTERNAL FP 

1F(NT.NE.30)GO TO 50 

R2=(R-1.0D0)«i«2 

KC=1 

2 A2=i.0D0-B 
XZ=RM**A2 

AU=AZ4D5QRT( (X-XZ )«^»»24RZ ) 

1F(X.LT.XZ)G0 TO 10 
S=KOOO 
5 AL=AZ 

ANM=ZER0(FP,AL,AU,1 ,00-5, EPS) 

AMN=AMNf D5QRT( ( X-RM*<AMN-S»»DSQRT( AMN>»AMN-AZ*AZ ) ) 

GO TO 20 

to IFCRM.GT.1 .ODO)GO TO 15 
G=-t.0D0 
GO TO 5 
15 AMN=AU 

20 IF(NC.EQ.2)GO TO 22 
T(2)=AMfl 
T(1)=AMN-1.0D-^ 

T(3)=AttrUK0D-^ 

B=-B 
HC=2 
GO TO 2 
22 T(5)=AMM 

T(<<)=AMN-1 .00-^ 

T(6)=AMNfl .00-^ 

T(7)=DINT( At:N+l .ODD) 

DO 25 1=6, n 
25 T(I)=TtI-l l^l.ODO 
DO 30 I=Kf,23 
30 TCI)=T(r-l )+2.0D0 
TE=T(23) 

EL=30.0DC 
DO 35 J=T,20 
1F(TE.LT.EL)G0 TO ^0 
35 EL=EL45.0D0 
^0 TC2^^ = EL 

DO ^5 1=25.29 
A5 T(I)=T(I-1 )45.0D0 
T(30)=T(29)+10.0D0 
D=~B 
RETURN 

50 READ 55,(T(I),I=l ,NT) 

55 FOPMAT(7D10.3) 

RETURN 

END 

C I* « X M K M ¥ X M « K » » |i! M K K » K M M ^ M ^ H )» M ^ .1 

FUNCTION FP(T) 

C 

C FUNCTION EXTEPlJAL TO TIME AMO CALLED BY ZERO 
C 

IMPLICIT PEAL*‘a(A-»l,0-2> 
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COMMON /CMA/ RM,R,B»A,X /CTA/ AZ^5»R2 
AArD5QRT(Ti»T-A2«^A2) 

AF=X-RM»»T-5»fAA 

FP=DSQRT( AF*iAF4RZ)<<AA-AF«<(RM^AA*S«T) 

RETURN 

END 

SUBROUTiriE ROOT(U,V,ISP,FS) 

c 

C PROGRAM FINDS THE ROOT OF F OR G 
C FA,ZA “ GAMMA ONE. FDiZB - GA?1MA TWO 

C FS - VALUE OF F OR G AT THE ROOT 

C 

IMPLICIT REAL*f0(A,B,D-H,0-Y)*C0MPLEX»^16(Z) 

COMPLEX^t6 CDSQFU.CMPLX 

DIMENSION BrH5),FR(^),FI(A),FMA(^),FMB(^),FMS(^),XM(^) 
CO^:MO^i /CMA/ RM,R,B,A,X /CRA/ FA.FB /CRC/ JS /CRB/ ZA*ZB 
JS=0 

RMS=RM^RM 

D5={ 1 .0D0-D)»^>^2 

AP=2.OD0^»(RM‘^B5-X*fA) 

AZ=X*^X4B5^( 1 .0D0-RM5) 

AT=A 4 fA-B 5 

IF(R.EQ. 1 .ODO)GO TO 20 

R5=(R-1 .0D0)K^(2 

AS=A.0D0»*R5«B5 

AZ=AZ4RS 

AT=AT-RS 

BM(5)=AT*<AT-AS 

BriCA) = (AP>‘AT4A54<RM)4f2.0D0 

BM(3)“AP«AP+2.0DO^‘AZ**AT-(RM5-2.0DO )»*AS 

DMf 2 ) = 2 , 0D0i4 ( APHA2-AS^<RM ) 

BMt 1 )=AZ4^AZ4AS<4(RMS-1 .ODO) 

CALL QUARTKBM.FR.FI) 

JD=A 
5 MM=0 

DO 15 J=t ,JD 

1F(FI( J).EQ.O.ODO)GO TO 10 
IF(FI(J).LT.O.ODO)GO TO 15 
2=DCMPLX(FR(J),F1(J)) 

ZA=CDSQRT(Z4fZ“l .ODO) 

ZB=CDSQRT( (Z-RM)K4‘2-1 .ODO) 

IF ( DIMAGt ZA ) . LT. 0 . ODO )2A=-ZA 
IF ( DIMAG( ZD ) . LT. 0 . ODO )2D~-ZB 
ZS=X4ZA^(R-1 .0D0)4ZD»4( 1 .ODO-Bl-A^'Z 
FS=CDABS(ZS) 

IFCFS.GT.I .OD-A)GO TO 15 
1SP = 2 
U=FR( J) 

V=FI( J) 

RETURN 

e FA=DSQRT( 1 .0D0-U<»U) 

ZA=DCMPLX(O.ODO,FA) 

2B=0CNPLX(FD,0.0D0) 
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JS=1 

CO 70 12 
10 U=FR(J) 

FB-DSQRT( (U-Rm»»»»2-1 .000 ) 

1F(U.LT.RM-1 ,0D0)FB=-FB 
IF(DABS(U).LT.1.0D0)GO TO 6 
rA=DSQRT(U»U-t .000) 

IF(U.LT.-1.0D0)FA=-FA 
12 FS=X+FA>^(R-1 .ODO) + FB'(( 1 .0D0-B)-A^U 
IF(DADS(F5).GT.1.0D-^)GO TO 15 
m=m*i 
XMt«M)=U 
FnA(Mm=FA 
FM3(riM)-FB 
FHS(m)=FS 
15 CONTINUE 

CALL SELECT(m,XN,MN) 

FA=FMA(NN) 

FD=FnB(NN) 

U-XM(NN) 

V=O.ODO 

FS=FnSfMN) 

IF(ISP,EQ.2)ISP=3 

RETURN 

20 AA=AP*‘AP-A.0D0>^AT*^A2 
SAD=-AP/AT^0,5D0 
IF(AA.LT.O,ODO)GO TO 25 
saa=D5»:;rt(aa)/at^o.5Do 
FRM )=SAD + 5AA 
FR(2)=SAD-SAA 
FI(t ) = 0.0D0 
FI(2)=0.000 
JD=2 
GO TO 5 

25 5AA=DSQRT(-AA)/AT^*0.5D0 
FR( 1 )=SAD 
FR(2)=5A0 
FI( 1 )=SAA 
FI(2)=“5AA 
JD=2 
GO TO 5 
END 

FUNCTION 5PEC(U»V) 

C 

C CALCULATES THE SPECULARLY TRANSMITTED FIELD 

C U AND V ARE RESPECTIVELY THE REAL AND IMAGINARY PART OF PSI 

C 

IMPLICIT REAL»^G(A,B,D-U,0-Y), C0MPLEXxi6( C ,Z ) 

COMMON /CMA/ RM,R,B,A,X /Cf.B/ DUM.FP.FI /CMC/ NN /CRA/ FA»FB 
COMMON /CRC/ J5 /CRD/ ZA,ZD 
RR=R-1 .000 
DD=1 .000-0 

IFCV.EQ.O.ODO.AMD.JS.EQ.O)GO TO 10 
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ZFAC=OCnPLX(FR,FI) 

Z=DCnPLX(U,V) 

ZG=(ZA»ZB)«f(ZA‘»ZO + l .ODD) 

ZM=Z*^(Z-RM) 

ZC=(Z-RM)^BD/ZB 

IF ( . HE , 3 )ZC=ZC*Z^«RR/ZA 

ZFH=ZH/ZG 

ZFD=ZC-A 

IF(HW.EQ.2)GO TO 5 
IFnW.EQ,3)2FN=-ZFH^<ZA/Z 
SPEC=DREAU ZFH/ZFDJtZFAC ) 

RETURN 
5 HW=1 

ZMI=DCMPLX( 0.000,-1 .000) 

ZFG=X*ZA*fRR>ZB»*BB 

ZDT=Z>^Z/(Z^^ZC-ZFG) 

ZDG=2.0D0*‘(Z^ZB*(2-Rrt)i<2A)+2rt«( (2-RM)/ZA*2/ZB) 
ZDQ=(ZGH(2.000*Z-Rri)-2M*fZDG)/(2G^^ZG) 

ZDFD = (ZA-2^2/ZA)»^RR/(ZAitZA)-Z0T 
♦ 4(ZD-tZ-RH)J^(Z-RM)/ZB)vDD/(ZB^^ZB) 

SPEC=-DREAL( tZFD>*Z0P-ZFH^ZDFD)/(ZFD^2FD)v2rAC-Zm»^FUNCB(Vn»ZDT 

uu-z 

RETURN 

10 FG=(FA*FB)»^(FA^FB+1 .000) 

FN=U^t(U-RM) 

FC=(U-RM)»«BB/rB 

IF(HW‘NE.3)FC=FC*U»^RR/FA 

FN=-FIt/rG 

fd=fc-a 

IF(NM.EQ.2)G0 TO 15 
IFaW.EQ.3)FN=-FN»*FA/U 
SPEC=FK/FD 
RETURN 
15 HW=1 

FFG=X+FA»‘RR+FB>tB3 

FDT=U»^U/(U^‘FC-FFG) 

FDG=2.000v*(U»‘FB*(U-RM)#FA) + FH*fC(U-Rri)/FA+U/FB) 
FDQ=(FG^(2.0D0>iU-RM)-FHitFDG)/(FG^FG) 

FDFD = ( FA-UJ^U/FA )^*RR/( FA^FA )-FDT 
4 4(FD-lU-Rn)*^(U-RH)/FB)ifBB/(FB^FB) 

SPEC=-FDTH( ( FDCJ4FD-FDFD*^FN )/( FD>JFD )-FUNCA(U) ) 

HW=2 

RETURN 

END 

FUNCTION SMOD(U) 

C 

C CALCULATES THE FIELD UMICH MODIFIES THAT OF THE 
C SPECULARLY TRAt^SMITTEO WAVE 

C P?ESE»iT BEHIND THE NEUTRAL STABILITY WAVE 

C 

IMPLICIT PEAL*‘0(A-H,0-Z) 

COMMON /CMA/ RH,R,D,A,X /CMC/ »W 
DATA T5T/2.O280C7t2^t7^619D0/ 
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SnOD=O.ODO 

inRM.LT.2,0D0)RETURN 

IFJ A.LT.RM»(R-B)/DSQRT(RM»RM-^.ODO))RETURH 
GA=DSQRT(U^U-1.0D0) 

GB=-GA 

PF=X+GA»*CR-1 .ODO)+GB^( 1.0D0-D)-A-U 
IF(RM.GT.TST) GO TO 10 
IF(PF.GT.O.ODO) RETURM 
FAC=1 .ODO 
5 GA3=GA«GD 
GPT=1 .ODOfGAQ 

QT=U*^(U-RM)/(GPT^^(U/GA+(U-RM)/GB)»<O.ODO) 

IF(NW.EQ.3)QT=-QT*tGA/U 

PS=tGA-G3)»^( 1 .0D0-3.0D0>^GAQ)/GPT/GAB 

PSM=PS“GA/D+GB/R 

IF(K’Vi.EC?.3)PSM=PS-GA/B-3,0D0**GD 

SMOD = P3M*fQT^‘DEXP(-P5*^PF/6.0D0)^^FAC 

RETURN 

10 IF(PF.LT.O.ODO)RETURM 
FAC=-1 .000 
GO TO 5 
END 

FUNCTION LAG(FUHC,IC,I,SP,MI,DI,ISP) 

C 

C PERFOPMS INTEGRATION OF FUNC OVER (5P,INF) 

C BY GAUSSIAN irUEGRATIONS. IF SP < 0 , THEN I»4TECPATES 

c OVER (SP,-irir) 

C RANGE IS SPLIT INTO HI PARTS. (SP»SP1), (SPWSP2), • 


C (SP(NI-1 ),ir^F) 

C THE riUMBER OF GAUSS POINTS RANGE FROM 

C MINinUfl OF 2 TO MAXIHUM OF 16 

C SP( J)-SP( J-1 )*DI( J) 

C SUPPLIED BY USER TO SUIT HEED 

C Y - GAUSS POINTS ARRAY LESS ALL ZEROS 

C V - GAUSS WEIGHTS ARRAY WITH THE ADDITION OF 

C VZ - GAUSS WEIGHTS ARRAY WHEN CORRESPONDING POINTS ARE O.ODO 

C JRG - POINTER TO Y AND V 

c jz - PoirncR to vz 

c 


IMPLICIT REAL*BtA-H»L,0-Z) 

DINENSICN JRGM2)»JZ( n ) i IDC5 ) >V( 39 ) , Yt 39 ) ,VZ( ^ ) ,DD( BO ) ,E( 5 ) ,DI( 1 ) 
DATA JRG/1 ,2,3,5,7,10,13»17,21 ,26>32,A0/ 

DATA JZ/Sfl ,5,2,5,3,5,^,5,5,5/ 

DATA Y/0.5773S0269189626t)0,0.77^5966692AlA03DO, 0.661 13631 159A053D0 

♦ »0. 3399S I 0^356*^65600,0. 9061 793^593C66AD0,0.538A67310105663D0 

♦ •0.932^*6951620315200,0.66120933660626500,0.23061918606319700 

♦ . 0.9691 0791 236275900,0. 761 5311 C559939;D0, 0.6058651 51 377397D0 

4 •0.96020965669753600,0.79666667761362700,0.52553260991632900 

♦ .0.18363666269565000,0.96016023950762600,0.83603110732663600 

♦ .0.61 337163270059000, 0.32625362360330900,0.97390652651717200 

4 •0.86506336660093500,0.67960956629902600,0.63339539612926700 

4 • 0.16007633896163100,0.98156063626671900, 0.90611 7256370675D0 

4 ,0.76990267619630500,0.50731795626661700,0.36783169099818000 
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♦ »0, 125233^06511^6900,0. 989^0093^991650DO,0.9'*<457502307323300 

4 , 0,665631 202367O32D0,O.755<»C^«»0e355003D0,0. 61 76762^*^^026'*<iD0 

4 ,0,^5001677765722700,0.2816035507792600,0.09501250933763700/ 

DATA V/l.ODO ,0.555555555555556DO,0.3^785^S'i5137^5^00 

4 ,0.652U515^0625<<6DO, 0.23692688505618900,0. <*78628670^9936600 

4 ,0.1 71 326<#923791 7000,0. 3607615730<iB1 3900, 0.<i6791393<*572691D0 

4 ,0.12948696616887000,0.27970539168927700,0.38183005050511900 

4 ,0.10122853629037600,0.22238103665337600,0.31370666587780700 

4 ,0.36268378337036200,0.08127638836157600,0.16066016069685700 

4 ,0.26061069660293500,0.31236707706000300,0.06667136630868800 

4 ,0.16965136915050100,0.21900636251598200,0.26926671930999600 

4 ,0.29552622671675300,0.06717533638651200,0.10693932599531800 

4 ,0.16007832656336600,0.20316762672306600,0.23369253653835500 

4 ,0.26916706561360300,0.02715265961175600,0.06225352393666000 

4 ,0.09515851166269300,0.12662397125553600,0.16959593801657700 

4 ,0.16915651939500300,0.16260361506692600,0.1696506106550700/ 

DATA V2/0.eOOSD86SB8806S9DO, 0.56860883608838900, 

4 ,0.61795918367366900,0.33023935500126000/ 

SGM=1 .000 

iniGP.EQ.l )SGN=D5IGN( 1 .ODO.SP) 

5V=SP 

1C=0 

DO 5 1=1, NI 
5 DD(I)=Dia) 

KL=1 
KLS=1 
1S=1 
IL=UI 
GSD=O.DDO 
10 JS=1 
JP=0 
JML=1 

DO 55 1 = 15, IL 
A=SP 

B=SP*DDm*<5Gtl 

q=0.5D0*^(D*A) 

P=Q-A 

TD=O.ODO 

IB=0 

DO 25 J=JS,11 

CD=O.ODO 

JF=JPGU) 

JL=JPG( J+1 )-l 
JM=J2(J) 

DO 15 K=JF,JL 

ID=IB*2 

Pr=p4Y(K) 

15 GD=G0W(K)*^(FUT4C(Q-Pr)*FUNC(Q*PY) ) 

1F( JM.EQ.5JG0 TO 20 
GD=CDfV2( Jm**FUNC(Q) 
lB=lDtl 
20 TA=CD»^P 

lF(DADSnTA-TB)/TA).LT.0.5D-3)GO TO 35 
25 TB=TA 
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PRINT 30»1>TA 

30 rORMATC •»*riAX*»12»WAL*,010.3) 

ID(KL)=I 
E(KU=A 
KL=KL*1 
GO TO ^2 
55 GSA=GSD*TA 
1C=IC*IB 

1F(I.LT.15)G0 TO ^0 

1F(DAB5((GSA-G5B)/GSA).LT.2.0D-3)GO TO 60 
AO GSD^GSA 
A 2 JT=J-jriL 

IF( JT)A5,A7,50 
A5 J5=J-3 

IF( JT.LT.-l )JS=J-2 
GO TO 52 
A7 J5=J-3 

IF(JT,EQ.JP)J5=J-2 
GO TO 52 
50 JS=JtJT-2 

IF( JT.GT.2)JS=J+1 
52 JP=JT 
JML=J 

IF( J5.LT.1 )J5=1 
IF( J5,GT.10)JS=10 
55 SP=B 
60 COtUIN’OE 

IF(KLS.EQ.KL)GO TO 70 

1P-ID(KLSJ 

SP=E(KL5) 

D5=DD(IP)/1 .001 
1S=1 
IL=I49 
KL5=KL5+1 
DO 65 JJ=IS,IL 
65 DD(JJ)=DG 
GO TO 10 
70 LAG=G5A 
SP-SV 
RETURN 
END 

M M K K X K K K M ^ ^ ^ ^ ^ 

FUNCTION rUNCAtU) 

c 

C INTEGRAL ALONG THE BRANCH CUT DEEDRE WAVE ARRIVAL 
C U - ARGUMENT SUPPLIED BY FUNCTION LAG, REAL PART OF PSI 

C 

ItlPLICIT REALx6(A-H,0-Z) 

COKHDN /CtlA/ RM,R,D,A,X /CMC/ NW 
PF=1 .ODD 

GA-DSttRKU'U-l .ODD) 

CD-DSORT( (U-RH .ODD ) 

IF(U.LT.-1 .0D0)GA--GA 
IF(U.LT.RM-1 ,0D0)GB=-GB 
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CACsGAfCB 

FP=n .ODO«GAD>«GAB 

SP=tGA-CB)x( 1 .ODO-3.0DO«GAB)/FP 

IF(MW.EQ.E)PF=UxSP/6.0D0 

TP=UX(U-RM)/(GA*GB) 

SD=TPxSP/t 1 .ODO+GAB) 
lF(N«.Eq.t )TO=TP»tGB/R-GA/B)/FP 
IF(liW.EQ.3)T0=TPxGA/U»(3.0D0xGD + GA/B)/FP 
AM=DI‘x(XtGBx( 1 .ODO-B)-A*U)/8.0DO 
IFtNM.NE.3)AM=AM+SPxGAx(R-l .ODO )/8.0D0 
FU«CA=-( SD+TD )XDEXP( AM )XPF 
RETURU 
EtJD 

FUUCTIOM FUS’CD(V) 

C 

C IIITEGRAL ALOMG THE BRAMCH CUT AFTER WAVE ARRIVAL 
C V - ARGUrlEHT SUPPLIED BY FUHCTIOH LAG, IMAGINARY PART OF PSI 

C USP- UO, STARTING POINT, REAL PART 

C FCR- 0 FOR PHIA, -1 FOR PHIB 

C FCI- -1 FOR PHIA, 0 FOR PHIB 

C 

IMPLICIT REALX8(A,B,D-H,0-Y), C0MPLEXxl6(C,2) 

COMMON /CMA/ RM,R,B,A,X /CMB/ USP, FCR, FCI /CMC/ NW 
2PF=DCMPLX( 1.000,0.000) 

2FAC=DCMPLXtFCR,FCI) 

Z=DCMPLXtUSP,V) 

IF(V.Eq.O.ODO)GO TO 10 
2A=CDSqRT(ZxZ-l .000) 

2D=CDGQRT( (2-KM)xx2-i .000 ) 

IF ( DIMAG( ZA ) . LT. 0 . 000 )ZA=-ZA 
IF( DlMAGt ZB ) . LT. 0 . ODO )ZB=-ZB 
5 ZAB=ZAXZD 

ZFP=( 1 .0D0+ZAB)xZAB 
ZSP=(ZA-28)x( I .000-3. 0D0XZAB)/2FP 
IF ( t!W. Eq . 2 )ZPF=ZX2SP/8 . ODO 
ZTP=2 x12-RM)/(ZA+ZB) 

25D=ZTPxZ5P/( 1 .ODOtZAB) 

IF( N;J. NE . 3 )ZTD=ZTPX( ZA/B-ZB/R )/ZFP 
IF(NM.Eq.3)ZTD=-2TPxZA/Zx(3.0D0x2B+ZA/B)/ZFP 
ZAM=ZSPX(X4ZBX( 1 .0D0-B)-AXZ)/6.0D0 
IF ( N’.J . HE . 3 )ZAM=ZAM4Z5PXZ AX ( R- 1 . ODO )/3 . ODO 
FUNCB=DREAL( ZFACx( ZSD-ZTD )xCDEXP( ZAM ) xZPF ) 

RETURN 

10 FD=D5QRT((USP-RM)x»2-1.0DO) 

IF1U.lt. RM-1 .0D0)FB=-FB 
2B=DCMPLX(FB, 0.000) 

IF(DABStUSP).LT.1 .ODO)GO TO 15 
FA=DSQRT(U5P»>U5P-1 .000) 

IF(U.LT.-1.0D0)FA=-FA 
2A=DCMPLX( FA, 0.000) 

GO TO 5 

15 FA-DSqRTt 1 .ODO-USPxUSP) 

2A:DCMPLX(0.0D0,FA) 
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GO TO 5 
END 

C « KK MM If « M If « ¥ K » M »r M X )f K K M # X K M if K K K 

5UBR0UTIME SELECK M*>CM,K ) 

C 

C PROGRAM FINDS THE ROOT EXTERIOR OF CFO OR CGO 

C 

IMPLICIT REAL«8(A-H,0-2> 

DIMEfiSIOW XM(^*),XM(^) 

COMMON /CRC/ JS 
IFtM.EQ.I )GO TO 30 
DO 5 1=1 ,M 
5 >04(I) = DABS(XM(in 

GO TO (30,10^15,20), M 
10 XN(3)=O.ODO 
15 XN£0)=O.ODO 

20 X=DMAX1 ( XN( 1 ) »XN( 2 ) ,XN( 3 ) ,XN( ) ) 

IFUS.EQ.I )X=DMINnXM(l ),XN(2)) 

DO 25 K=1 

1F(>:.EQ.XM(K))RETURN 
25 CONTINUE 
30 K=1 
RETURN 
END 

FUNCTION 2En0(F,AX,BX,T0L»CrS) 

C 

C PROGRAM FINDS THE ZERO OF A FUNCTION F 

C 

IMPLICIT REAL^^etA-lUO-Z) 

A = AX 
B = DX 
FA = F(A) 

FB = F(B) 

lF(D5IGNM.CD0,FA).NE.DSIGtJ(1.0D0,FB))GO TO 20 
ZERO=A 
RETURN 
20 C=A 
FC = FA 
D = B - A 
E = D 

30 IF (DADSirC) .GE. DABS(FB)) GO TO ^0 
A = D 
B = C 
C A 
FA = FD 
FD = FC 
FC = FA 

<*0 TOU = 2.0DO»^EPS^DABS(D) ♦ O.SDOxlOL 
XM=0.500 m(C-B) 

IF (DADStXM) -LE. TOL1) CO TO 90 
IF (FD .EQ. O.OUO) GO TO 90 
IF (DABS(E) .LT. TOLD GO TO 70 
IF (DADS(FA) .LE, DAD5(FB)) GO TO 70 
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IF (A .ME. C) GO TO 50 
S = FD/FA 
P = 2.0D0>fXH»»S 
Q = 1.0D0 - 5 
CO TO 60 
50 Q = FA/FC 
R = FB/FC 
5 = FD/FA 

P = S»(2.0D0^XM»tQ»(Q - R) - (B - A)^tn - t.ODO)) 

Q = CQ - 1.0D0)*<(R - 1.0D0)H(5 - t.ODO) 

60 IF (P .GT. O.ODO) Q = -Q 
P = DABS(P) 

IF ((2.0D0i*P) .GE. (3.0D0^»XM<fQ - DAB5( T0L1<«Q ) ) ) GO TO 70 
IF (P .GE. DAB5(0.5D0*^E>^Q)) GO TO 70 
E = D 
D = P/Q 
GO TO 60 
70 D = XM 
E = D 
60 A = B 
FA = FB 

IF (DABS(D) .GT. TOLD B = B 4 D 

IF (DABD(D) .LE. TOLD B = B 4 DSIGMCTOLI, Xrt) 

FB = F(B) 

IF ((FD*(FC/DAD5(rC))) .GT. O.ODO) GO TO 20 
GO TO 30 
90 ZERO=B 
RETURN 
END 

SUBROUTINE QUARTI (C,XR,XI) 

C 

C PROGRAM FI)iD$ ROOTS OF QUARTIC EQUATION 
C CD ) 4 CiZ)^X 4 C(3)^X*t*f2 4 C(^)>^X«it3 4 = 0 

C REAL PARTS OF ROOTS ARE PLACED IM XR, IMAGDiARY PARTS IN XI. 

C 

IMPLICIT REALMS (A-H,0-2) 

DinEM5IO>4 CCS), XR(A), XI(^), CU(5), YR(A), Yl(4) 

2=CC5) 

IF(2.KE,O.ODO)GO TO 20 
CALL CU3IC(C,XR,XI) 

XR(^)=O.ODO 

XI(^)=0.0D0 

RETURN 

20 XIM )=O.ODO 
XI(3)=0.0D0 
CU(5)=O.ODO 
CUC^») = 1 .000 
CU(3)=-C(3)/2 
A0=C(D/2 
A12=CC2)/Z«0.5D0 
A32=C(^)/Z<<0.5D0 
cue 2 )= ( A32#A12-A0 )^*A.ODO 

cue 1 )= -(A12><«24A0*»( A32«*<24CU( 3)) )*»^.0D0 
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CAU CUCIC(CU,YR,YI) 

DO <i0 1=1.3 

U12= YR(I)*t0,5D0 

TI=UI2«>i2-A0 

ifcti.lt. o.odoigo to <»0 

T2=A32«»24UI2+U12«CU(3) 

IF (T2.GE. 0.000)00 TO 50 
<i0 CONFItCUE 

XIC 1 )=-)00.0D0 
RETURN 

50 TI=DSqRT(Tl ) 

T2=DG1GN1 DSORTI T2 ) , A32»<U1 2-A1 2 ) 

60 Cq=U12-Tt 

BQ2= (A32-T2)*t0.500 
DlSC=B02»t»i2-CQ ' 

IF CDI5C) 90 , eo. 70 
70 DI5C=DSQRT(DISC) 

60 XR( 1 ) = -BQ2<DISC 
XRC2)=-BQ2-DISC 
GO TO 100 
90 XR1 1 )=-BQ2 
XR(2)=-B02 
Xin )=DSQRTt-DlSC) 

100 Cq=U12+Tl 
BQ2=BQ2+T2 
DISC=BQ2>^*<2-CQ 
IF (DISC) 130, 120, 110 
110 DISC^DSIRTCDISC) 

120 XR(3)=-Bq2+DISC 
XR(^) = -BP2-DISC: 

GO TO 150 
130 XR(3)=-BQ2 
XR(''.) = -BQ2 
XI(3)=D5PRT(-DISC) 

150 XI(2)=-XI(1 ) 

XIC^)=-XI(3) 

RETURN 

END 

SUBROUTINE CU31CC C,XR .XI ) 

C 

C PROGRAM FINDS ROOTS OF CUBIC EQUATION 

C C( 1 ) + CC2)X>: ♦ C(3)**XXX2 + C(<«)xXxx3 r o 

C REAL PART OF ROOTS ARE PLACED IN XR, IMAGINARY PART IN XI. 

C 

IMPLICIT REALMS (A-H.0-2) 

DIMENSION C(5),XR(A),XI(<t) 

A=3.0D0 

2=C(A)x3.0D0 

irCZ.NE.O.ODO )GO TO 20 

XR(3)=0.0D0 

X1(3)=0.0D0 

C3= C(3)x2.0D0 

CC=C(2)xx2- C( 1 )«C(3)xA.0D0 
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IFICC)10,S0>90 

10 5C=DSQRT(-CC) 

XI 1 1 J=SC/C3 
XI(2)=-XI( 1 ) 

XR(I )=-C(2)/C3 
XR(2)=XR(1 ) 

RETURN 

50 XR( I ) = -CI2)/C3 
XR(2)=0.0D0 
no XI( 1 )=0.0D0 
XI(2)=0.0D0 
RETURN 

90 XRM ) = (-C(2)+DSQRT(CC))/C3 
XR( 2 )=( -Ct 2 )-DSQRT( CC ) )/C3 

CO TO no 

20 XI( 1 ) = 0.0D0 
XI(2 ) = 0.0D0 
P3=Ct3)/2 
Q3=C(2)/2 
PQ16=P3*iP3<<0.5D0 
R2=Ct 1 )»t .5D0/Z 
AmN3=P3'(*i2-q3 
B2:ANIN3>fP3-PQ16tR2 
DISC:B2**»2-AHIN3)H*3 
IP (DISC) 80> <tO> 30 

30 D1SC=DSQRT(DI5C) 

AO B2D=-B2*DISC 

SI=DSIGN( 1 .000,020) 
CAPA-SI»DAES(D2D)«v( t .000/3. ODO ) 

IF (CAPA.NE.O.ODO)GO TO 60 
CAPB=(-B2-B2)*»(1 .000/3.000) 

GO TO 70 

60 CAP3-AHIN3/CAPA 

70 XR( 1 )=CAPA+CAPB-P3 

XRt2) = -(CAPA+CAPB)**0.5D0-P3 
XR(3)=>:R(2) 

XIC2)=D5QRT(A ) «(CAPA-CAPB )tO .500 
CO TO 100 

60 R00T=D5QRTt AM1N3) 

PHI3=OARCOS(-B2/(ROOT*^AmH3))/3.0DO 
TnOSORTCA )*<D5IN(PHI3) 

T2=DC0S1PHI3) 

XRM )=T2>!ROOT>52.0DO-P3 
XR(2)=ROOT-(-T2-TI )-P3 
XR(3)=RC0T«t-T2+T1 )-P3 
100 XI(3)=-XI(2) 

RETURN 

END 
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